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VsiMJ^RY 

Research  efforts  covered  by  the  subject  contract  have  been  directed 
toward  the  stability  and  synthesis  problems  for  two  dimensional  digital 
recursive  digital  filters.  Such  filters  are  advantageous  when  a premium 
is  placed  upon  computer  storage  or  computation  time. 

Previously  defined  results  on  stability  analysis  have  been  refined  and  a 
revised  paper  has  been  submitted  to  the  IEEE  Transactions  on  Circuits  and  Systems 
for  publication.  This  research  effort  will  continue  to  address  this  complex 
problem. 

-'■'previously  defined  results  on  synthesis  of  bandpass  and  band  enhancement  filters 

/ 

/ have  been  refined.  Results  of  these  improvements  will  be  installed  in  the  Naval 
/ Intelligence  Support  Center's  Spatial  Domain  Filtering  Package  which  was 

designed  and  implemented  by  this  author.  This  work  will  be  done  under  a separate 


refinements  should  significantly  improve  the  overall  performance  of  the 
Documentation  on  this  spatial  domain  filtering  package  should  be 


I l effort.  This  package  has  provided  very  satisfactory  results  to  date  and 
| V the  refi 
\ package. 

| available  within  the  next  year. 

Approximately  circularly  symnetric  lowpass,  highpass,  bandpass,  band- 

I 
I 
I 
I 
l 
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stop,  low  frequency  boost  and  high  frequency  boost  filters  have  been  designed. 
Evaluation  of  these  filters  with  actual  images  of  various  types  need  to  be  addressed 
and  has  been  proposed  as  a follow  on  to  the  effort  described  in  this  report. 


l-wmlj-  ID 

re  me 

m Ma 

IIUMMCO 

jBTiFicinaa 

r* 

□ 

|*pe<  VW. 

. 

= 

| NTHNttM/lfJIUIIUTT  M8  | 

ML 

P 

3 


INTRODUCTION 

The  two-d iroens ional  recursive  digital  filter  is  particularly  suited  to 
image  processing  applications  when  there  is  a premiun  placed  on  computer 
memory  requirements  and  time  for  processing.  Due  to  these  considerations,  it 
has  definite  advantages  over  the  Fast  Fourier  Transform  (FFT)  algorithm  for 
many  image  processing  operations [l].  The  application  of  recursive  digital  filters  to 
image  processing,  however,  has  been  hampered  by  two  problems:  stability  and 
synthesis  [2 ]•  The  synthesis  problem  is  the  problem  of  expressing  the  desired 
impulse  response  in  closed  form  and  thus  detenu ing  the  filtering  coefficients. 

The  stability  problem  occurs  because  the  recursive  filter  requires  feedback 
of  past  output  values  and  therefore  it  can  become  unstable. 

The  research  for  the  first  year  has  been  on  both  problems  with  progress  made 
in  both  areas.  This  report  discusses  the  progress  made  in  both  areas  and  the 
directions  for  future  research. 

STABILITY 

The  stability  problem  for  one  dimensional  digital  recursive  filters  is 
straight  forward.  The  roots  of  the  polynomial  in  the  closed  form  of  the  one 
dimensional  2-Transform  for  the  filter  inpulse  response  must  have  magnitudes 
less  than  one.  Stability  analysis  is  therefore  reduced  to  finding  roots  of 
nth  degree  polynomials  with  constant  coefficients.  For  the  two  dimensional 
problem,  stability  is  not  straightforward  because  a two  variable  poly- 
nomial is  not  generally  factorable  into  distinct  roots.  When  the  polynomial 
in  the  denominator  of  the  two  dimensional  Z-Transform[3]  for  the  impulse 
response  is  factorable,  then  the  stability  analysis  procedure  is  the  same 
as  for  the  one  dimensional  problem. 

The  tvro  dimensional  stability  problem  is  very  complicated  if  the  poly- 


4- 


I 

l 

I 


T 

«*• 


I 


r 

1* 


nomial  in  the  denominator  of  the  Z -Transform  is  not  factorable  into  distinct 
roots.  Efforts  by  other  researchers  have  been  directed  toward  examining 
regions  of  roots  for  two  variable  polynomials  which  is  computationally 
feasible  only  for  very  simple  filters [4], 

The  method  used  by  this  researcher  is  to  express  the  two  dimensional 
digital  recursive  equation  as  a matrix  recursive  equation.  The  description 
of  the  matrix  recursive  equation  and  its  derivation  is  given  in  Appendix  C. 

The  resulting  matrix  recursive  equation  has  three  coefficient  matrices,  B^, 

B2  and  A.  Appendix  A gives  a surmary  of  stability  analysis  results  to 
date[5].  A paper  entitled  "Stability  Analysis  of  two  dimensioned  Recursive 
Filters"  by  W.  E.  Alexander  and  s.  A.  Pruess  was  revised  as  a part  of  this 
research  effort  and  resubmitted  for  publication  to  the  IEEE  Transactions  on 
Circuits  and  Systems.  A preprint  of  this  paper  is  given  in  Appendix  D. 

In  practice,  the  stability  analysis  procedure  which  only  involves 
finding  the  spectral  radius  of  a matrix  with  real  coefficients  is  very 
single  and  easily  implemented.  Computer  algorithms  are  readily 
available  to  perform  the  necessary  computations.  The  procedure  is 
regularly  used  by  this  researcher  for  stability  analysis  of  two  dimensional 
recursive  digital  filters. 

SYNTHESIS 

Often  it  is  possible  to  express  a desired  two  dimensional  digital 
recursive  filter  as  the  product  or  sun  of  two  one  dimensioned  digital 
filters.  That  is  the  two  dimensional  Z -Transform  of  the  digital  recursive 
filter  can  be  expressed  as  the  product  or  sun  of  two  one-dimensional  Z -Transforms. 
In  either  case,  the  tve  dimensional  synthesis  problem  is  reduced  to  the  synthesis 
of  two  one-dimensional  filters.  However,  it  is  not  possible  to  design  sun 
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separable  or  product  separable  digital  recursive  filters  for  all  applications. 
For  those  applications,  the  design  of  the  required  two  dimensional  digital 
recursive  filter  is  considerably  more  complicated. 

Many  imaging  systems  have  a natural  circular  synmetry.  In  general, 
the  optical  transfer  function  of  a circularly  synmetric  imaging  system  is 
uniform  with  respect  to  direction.  The  natural  consequence  is  that  filters 
with  circularly  symmetric  impulse  response  functions  are  generally  very 
desirable  for  image  processing.  The  relationship  between  circular  symmetry 
of  the  impulse  response  and  the  frequency  response  dictates  that  the  design 
requirement  is  for  these  filters  to  have  a circularly  symmetric  frequency 
response [6] . 

DOW  PASS  FILTER  DESIOJ 


The  design  goal  is  to  design  a low  pass  filter  with  circularly  symmetric 
frequency  magnitude  characteristics.  No  attempt  is  made  to  control  the  phase 
response  of  the  desired  filter.  This  presents  no  difficulties  in  im- 
plementing the  designed  filters  because  the  two  pass,  linear  phase  recursive 
digital  filtering  procedure  can  be  used  to  obtain  linear  phase [5] . 

The  magnitude  characteristic  for  the  one  dimensional  Butterworth 
approximation  filter  in  the  Laplace  Transform  variable  is  given  by 


h(s)h(-s) 
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The  corresponding  equation  for  two  dimensional  filters  is  given  by 
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where  s ^ and  <*>x  are  respectively  Laplace  Transform  and  cutoff  frequency 
variables  for  the  x direction  and  s 2 and  <*>y  are  respectively  Laplace 
Transform  and  cutoff  frequency  variables  for  the  y direction. 

If  the  bilinear  transformation [7]  is  applied  to  (2)  to  obtain  a two 
dimensional  Z -Trans form,  we  obtain 


H(z,v)2  - 


2 2 n 

[ (z+1)  (w+I)  ] 


[ (zH)2(w+i)2]n  + e2  (-l)nCn[  (z-1)  2 (w+1)  2+  (z+1)  2 (w-1)  2] 
' 1 /[tan2(wRT/2)]  u>RT  - ux2Tx2+uy 2Ty2 


Note  that  wy  is  the  effective  radial  cutoff  frequency.  In  continuing  the 
design  procedure  in  manner  similar  to  that  used  for  one  dimensional 
digital  recursive  filters [8]  difficulties  are  encountered  because  the 
denominator  of  (3)  is  not  factorable  in  distinct  roots  of  z and  w.  However, 
a suitable  approximation  may  be  obtained  by  factoring  along  the  w-z  plane. 


Thus  in  this  plane,  one  obtains 


H (z, z) 


(z+1) 


(z+1 ) 4n+ (-1 ) nCn [2 (z-1) 2 (z+1) 2] 


Simplifying,  we  obtain 


H (z,  z) 


(z+1) 

(z+1) 2n+e2 (-l)ncn[2  (z-1) 2]n 


Thus  the  poles  of  the  magnitude  response  in  the  w»z  plane  occur  in  reciprocal 
pairs  as  roots  of  the  denominator  of  (6) . 

As  with  one  dimensional  filters,  consideration  should  be  given  to 
round  off  errors  and  truncation  errors  in  implementing  two  dimensional 
digital  recursive  filters.  Thus  a cascade  realization  is  very  desirable 


because  it  acts  to  minimize  round  off  error.  Also,  it  is  desirable  to  avoid 
using  complex  arithmetic  when  implementing  two  dimensional  recursive  filters. 
This  leads  to  a natural  selection  of  implementing  a basic  filter  with  either 
one  pole  and  one  zero  or  two  poles  and  two  zeros  to  accomodate  complex 
conjugate  pairs  of  poles.  Then  any  general  filter  would  be  implemented  as 
stages  of  the  one  pole  or  two  pole  filter. 

If  we  let  n®l  in  (6),  we  obtain  a factorization  of  [h(z,z)|  2 * 
such  that  for  a stable  filter  design. 


H (z)»A  (z+1) 
(2+P) 


(7) 


As 


' 1 

I=3cF2 
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P - +(l+2Ce2  ) - 2^2C  e2 

1-2CE 2 (9) 

except  that  P*0  for  00.5/ e£ 

For  the  case  where  n is  equal  to  or  greater  than  2,  factorization 
becomes  more  complicated  and  the  computer  is  used  to  find  the  roots  with 
magnitudes  less  than  one. 

Forming  the  two  dimensional  Z -Transform  for  the  final  low  pass 
filter  design  for  n equal  to  one,  we  obtain 

H.(z,w)  - A2  (z+1) (wtl) 

(z+P)  (v+P)  (10) 


Note  that  this  filter  design  is  product  separable  and  inherently  stable 
because  we  have  selected  P such  that  |p|  is  always  less  than  one.  In  a 
similar  fashion,  we  can  design  filters  for  n greater  than  one. 

BAND  PASS  FILTER  DESIGN 

Once  a low  pass  filter  has  been  designed,  it  is  possible  to  obtain 
highpass,  band  pass  and  band  stop  filters  as  well  as  low  frequency  boost 
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and  high  frequency  boost  filters  from  the  low  pass  design.  In  this 
section , we  discuss  the  design  of  a general  boost  filter  which  can  be  used 
with  proper  parameter  values  to  obtain  the  above  mentioned  filters.  With 
the  low  pass  filter  design  (n»l)  given  in  (10) , we  can  obtain  a filter 
with  the  desired  magnitude  as  given  by 

| H(z,w)j  * a -*[hl(z,w)|2  (11) 

where  HL(z,w)  “HjUrW)  ^(z-3-,#"1).  Thus 
H(z,w)  * a + fi  Pi  (z+1)  (z-3*!)  (w+  1)  (w  +1) 


(z-FPI  fz-l+P)  (w+P)  (W"l+P)  (12) 

I H(z,w)l  »g[  (z+p)  (1+Pz)  (w+P)  (1+Pw)  1 +6rf  [(z+l)2(w+l)2] 

(z+P)  (1+Pz)  (w+P)  (1+Pw)  (13) 

Note  that  the  filter  represented  by  (13)  is  unstable  because  of  the 

terms  (1+Pz)  and  (1+Pw)  in  the  denominator.  The  poles  corresponding  to 

these  terms  have  magnitudes  greater  than  one  since  P has  a magnitude  less 

than  one.  However,  we  can  stabilize  (13)  by  using  the  minimum  phase  version 

of  the  denominator.  This  does  not  change  the  magnitude  since  the  magnitude 

of  (l4Pz)  is  equal  to  the  magnitude  of  (a+P)  and  the  magnitude  of  (1+BO  is 

equal  to  the  magnitude  of  (w*P) . Thus  the  desired  boost  filter  design 

has  the  two  dimensional  Z-Transform 

H (z,w)  - a [Pz2+  (1+P2)  z+P]  [P  w2+  ( 1+P2 ) w+P 3 + BA4r<z+l)2(w+l)21 


(z+P) 2 (w+P) 2 


thus  if  express  (14)  in  the  form 


Hfi(Z,W) 


z -Jw-K 


Z-JW“K 


it  follows  that 


gP2  + BA4 
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a02  “ a20  * ap2  + 8A4 
all  “ a (1+P2) 2 + 4 BA4 
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It  only  remains  to  determine  the  value  of  e for  both  the  low  pass  and 
boost  filters.  Note  that  the  squared  magnitude  of  Hl(Z,W)  is  equal  to 
1/(1+ when  the  radial  frequency, u , is  equal  to  the  radial  cutoff 
frequency,^.  Thus  for  n=l  |hl(Z,W^2  * l/d+e2)  at  the  cutoff  frequency. 

If  we  use  the  double  pass  linear  phase  filter  which  is  desirable [5] , the 
magnitude  of  the  resultant  filter  is  the  squared  magnitude  of  the 
orginial  filter.  If  we  desire  the  magnitude  of  the  resulting  filter  to  be 
down  3db  at  the  cutoff,  we  obtain 

1 » 1 

Y1  (l+e2)2  (18) 

The  same  value  for  e2is  appropriate  for  the  low  frequency  boost  filter. 

For  the  high  pass  filter,  we  have  a=l  and  fj*-l.  Thus  the  magnitude 
of  the  frequency  response  for  the  double  pass  linear  phase  filter  at  the 
cutoff  frequency  is  given  by  the  relationship 

1 - 1 - _1_  2 (19 

2 (1+e2) 
e2  - 1.0/(2*-l) 

The  same  value  for  £2is  appropriate  for  the  high  frequency  boost  filter. 

If  we  designated  B as  the  magnitude  of  the  desired  boost,  then  the 
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low  frequency  boost  filter  values  for  a and  8 are  given  by 

ot  = 1.0  ; 3 * B-1.0  (20) 

Correspondingly,  for  the  high  frequency  boost  filter,  the  values  of  a and 
0 are  given  by 

a = B ; 3 « - B+1.0  (21) 

Examples  of  two  dimensional  recursive  filter  designs  are  given  in 
Appendix  B. 

ROTATED  ONE  DIMENSIONAL  FILTERS 

A problem  of  interest  in  image  processing  is  to  filter  with  a one 
dimensional  filter  with  the  orientation  of  the  filter  specified  and 
•independent  of  the  sampling  directions.  This  type  of  filter  would  be 
useful  for  enhancing  or  suppressing  linear  features,  for  system 
noise  suppression  or  for  image  correction  (i.e.  linear  smear).  However, 
any  one  dimensional  digital  recursive  filter  which  is  rotated  becomes  a two 
dimensional  filter  associated  problems  in  stability  and  synthesis. 

Constraints  with  regard  to  angle  of  rotation  and  stability  of 
rotated  filters  have  been  developed  by  Costa  and  Ventsonopoulos[9] . They 
have  used  severed  rotated  low  pass  filters  to  obtain  circularly  symmetric  lowpass 
filters.  This  approach  is  currently  being  evaluated  with  regard  to  use  with 
single  rotated  filters. 

RECCMMENDATIONS  FOR  FITTORE  RESEARCH 

Approximately  circularly  symmetric  band  pass  and  band  enhancement 
filters  have  been  designed.  These  filter  designs  must  be  evaluated  with 
regard  to  performance  on  actual  images  of  various  types  and  with  regard 

H 
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to  circular  symmetry  as  a function  of  critical  frequency.  Methods  of 
improving  response  and  eliminating  errors  in  circular  symmetry  need  to  be 
investigated. 

Techniques  for  obtaining  rotated  one  dimensional  digital  filters  need 
to  be  investigated  with  regard  to  practical  use.  The  most  practical  method  will 
be  developed  into  an  algorithm  for  image  processing. 

The  use  of  recursive  digital  filters  for  image  correction  has  been 
hampered  by  problems  in  designing  filters  with  arbitrarily  specified 
magnitude  and  phase  characteristics.  Long  term  future  research  efforts  need  to 
be  directed  toward  this  problem. 

Recursive  digital  filters  are  practical  for  image  processing 
applications  using  small  computers,  minicomputers  and  special  signal 
processors.  Applications  range  from  real  time  image  processing  and 
data  acquisition  to  medical  applications  and  industrial  process 
monitoring.  However,  special  algorithms  must  be  developed  for  many 
of  these  applications.  Such  practical  research  problems  provide  the 
real  payoff  for  research  on  two  dimensional  digital  recursive  filters. 


APPENDIX  A 

Summary  of  Two  Dimensional  Digital  Recursive 
Digital  Filter  Stability  analysis  Results 


Given  the  two  dimens ionaldigital  recursive  filter  with  the  corresponding 


biavariate  recursive  equation 
L L 

g(m,n)  = Y1  'll 
J*0  K-0 


L L 


aJRf (m-J,n- 


n-K)  “ ^ bJKg(m-J,n-K) 

i=*0  K=0 


J+K>0 


where  g (m,n)  is  the  current  output,  g(m-J,n-K)  represents  past  output  values  and 
f(m-J,n-K)  represents  current  and  past  input  values  for  all  permissable  values  of 
J and  K.  It  possible  to  represent  this  relationship  by  a matrix  recursive 
equation  (see  Appendix  C) 

Gm,n  * BlGm-l,n  + *2Gm,n-l  + AFra,n 

where  G is  a column  vector  such  that  all  its  elements  are  the  outputs, 
m,n 

g (m-J,n-K) , where  0 and  0 L.  Fm>n  is  a column  vector  such  that  its 
elements  are  the  inputs,  f(m-J,n-K)  and  B j^and  A are  appropriate  coefficent 
matrices  such  that  (A.l)  and  A. 2)  are  equivalent. 

Theorem  1:  Given  the  discrete  system  represented  by  the  matrix  recursive 
equation  in  (A.2).  If  either  spectral  radii,  p(B^)  or  p(B^  is  greater  than  or 
equal  to  one,  then  the  system  is  computationally  unstable. 

Theorem  2:  Given  the  discrete  system  represented  by  the  matrix  recursive 
equation  in  (A.2).  The  system  is  computationally  unstable  if  p(Bi+B2)  *s  greater 
than  or  equal  to  one. 

Theorem  3:  Given  the  discrete  system  represented  by  the  matrix  recursive 
equation  in  (A.2).  The  system  is  stable  is 


£ abs  (Bj 


) + abs 


<i2)] 


(A. 3) 


where  abs(%)  and  abs(Bj)  refers  to  taking  the  absolute  value  of  all  elements 
in  and  B2  respectively  . 

Theorem  4:  Given  the  discrete  system  represented  by  the  matrix  recursive 
equation  in  (A.2).  Define  a particular  permutation  matrix  S (See  Appendix  C) . 
The  system  is  stable  if  P(Bi+B2)<1,  p(BiS)<l/2  and  p(%S)<l/2. 

Conjecture:  Given  the  discrete  system  represented  by  the  matrix  recursive 
equation  in  (A.2).  The  system  is  stable  if  and  only  if  p (B^)<1,P(B2)<1  and 
0(B1+B2)<1. 

Proofs  for  Theorems  1 through  4 have  been  developed.  Current  research 
is  directed  toward  verifying  the  practical  usage  of  these  theorems  and  to 


further  investigation  of  the  conjecture. 
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Filter  Design  Examples 

The  filter  synthesis  procedure  for  designing  two  dimensional  digital 
recursive  filters  in  this  research  effort  is  an  extension  of  a one  dimensional 
filter  synthesis  procedure.  The  squared  magnitude  characteristic  of  the 
desired  circularly  symnetric  two  dimensional  filter  is  chosen  in  the  Laplace 
Transform  domain.  The  Butterworth  filter  characteristic  has  been  chosen 
because  of  its  wide  spead  use  in  band  pass  filter  applications.  The  bilinear 
transformation  is  then  used  to  map  the  squared  magnitude  characteristic  into 
the  two  dimensional  Z-Ttansfocm  domain.  The  coefficient  matrices,  B 1 and 
B2  of  the  corresponding  matrix  recursive  equation  (See  Appendix  C)  are  obtained 
and  the  eigenvalues  of  the  matrix  sum  (B^  + B2)  are  determined.  These  eigen- 
values occur  in  reciprocal  pairs  because  the  original  function  was  a magni- 
tude response.  The  eigenvalues  with  magnitudes  less  than  one  are  then  used 
as  roots  of  a product  separable  denominator  to  form  the  denominator  of  the  two 
dimensional  Z-Transform  for  a stable  filter.  The  numerator  for  the  filter  is 
retained  from  the  mapping  of  the  squared  magnitude  characteristic  to  the 
two  dimensional  Z-Transform. 

This  procedure  has  been  used  to  design  and  implement  two  dimensional 
recursive  lowpass,  highpass,  low  frequency  boost  and  high  frequency  boost 
filters.  Some  examples  are  given  below.  It  should  be  emphasized  that  this 
design  procedure  always  results  in  stable  filters.  Also,  the  filter 
algorithm  necessary  to  implement  these  filters  has  been  developed  for  a CDC 
6400  System  by  this  researcher  at  the  Naval  Intelligence  Support  Center  in 
Suitland,  Maryland. 

There  are  still  some  minor  problems  remaining  with  this  procedure  with 
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regard  to  obtaining  circular  symmetry  when  the  critical  frequency  is  near  the 
Nyquist  frequency  or  near  zero.  These  problems  are  being  studied  with  the 
intent  of  providing  necessary  corrective  improvements  in  these  areas. 

Preliminary  results  indicate  that  the  problem  near  the  Nyquist  frequency  is  caused 
by  the  mapping  of  the  squared  magnitude  characteristic  to  the  two  dimensional 
Z -Transform  using  the  bilinear  transformation.  However,  it  does  appear  that 
significant  correction  can  be  made  for  the  case  when  the  critical  frequency 
is  near  zero. 

Figure  B.l  shows  the  contour  plot  of  a low  pass  filter  design  with  a 
cutoff  frequency  of  0.4  of  the  Nyquist  frequency.  The  contour  labeled  D is 

u 

the  half  power  point.  Figure  B.2  is  the  contour  plot  of  a high  frequency  boost 
filter  with  a break  frequency  of  0.5  and  relative  boost  of  high  frequency 
to  25.6.  The  contour  labeled  D is  the  half  power  point.  Figure  B.3  is 
The  prospective  plot  of  this  filter.  The  design  goal  is  that  the  contours 
and  specifically  the  break  frequency  contour  be  circularly  symmetric. 


20 


APPENDIX  C 

Description  of  the  Matrix  Recursive  Equation 
C.l  Introduction 

In  this  appendix,  the  matrix  representation  of  the  two  dimensional 
digital  recursive  filter  is  presented  in  detail.  The  matrix  S which  is  used 
for  the  Proof  of  Theorem  4 for  stability  analysis  as  presented  in  Appendix  A is  also 
described. 

C.2  The  Matrix  Recursive  Equation 

Consider  the  two  dimensional  digital  recursive  filter  which  has  the 
ZW-Transform  


I I 


-J  -K 
,z  w 


H(z,w)  - J~°  K“0 

L L 


(C.l) 


I z 

J-0  K-0 


, -J  -K 
b z w 
JK 


The  corresponding  two  dimensional  recursive  algorithm  for  the  filter  is 
given  by 


L L 


L L 


i,n)  - Z Z ajKf (m_J'n-K)  - Z Z bjK9(m-J,n-K) 


J-0  K-0 


J-0  K-0 

J+K  > 0 


(C.2) 


Define  the  matrix  V such  that  the  element  of  V in  the  Jth  row  and  Kth 

col um  is  given  by  b . That  is 

J“1,K-1 


b01 

• 

• 

• 

b0L 

bll 

• 

• 

• 

blL 

• 

• 

• 

• 

• 

bL0  bLl 


b. 
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or 


V “ [VJk1  “ [bJ-l,K-l] 


(C.4) 


V is  the  filtering  matrix  corresponding  to  the  denominator  of  the 
ZW-Transform  H(z,w) . A similar  filtering  matrix  can  be  defined  for 
the  numerator  of  the  ZW-Transform.  However,  this  matrix  is  not  used  in 
this  development. 

Define  the  vector  G which  contains  all  of  the  outputs  in  the 

m,n 

recursive  algorithm  for  the  filter  as  given  in  (C.2).  Order  the  outputs 

such  that  the  outputs  g(m-J,n) , for  all  values  of  J,  occur  before  the 

output  g(m,n-l)  and  the  outputs  g(m-J,n-l),  for  all  values  of  J,  occur 

in  order  before  the  output  g(m,n-2).  Continue  in  this  manner  until  all 

outputs  are  included.  The  vector  G is  then  given  by 

m,n 

g (m,n) 
g(m-l,n) 


G 


g(m-L,n) 

g(m,n-l) 


(C.5) 


9 


(ro-L,n-l) 


g(m-L,n-L) 


-J U.z 


- 


The  vector  G^  ^ is  then  obtained  by  decreasing  the  first  parameter  of 


each  output  in  by  one.  That  is 

g(m-l,n) 

g(m-2,n) 


m-lfn 


g(m-L,n) 

g(m-l,n-l) 


(C.6) 


lg(m-L-l,n-l) 


lg(m-L-l,n-L) 


Hie  vector  Gm>n_].  is  obtained  by  decreasing  the  second  parameter  of  each 


of  the  outputs  in  G . Thus 

m,n 


g(m,n-l) 

g(m-l,n-l) 


g (m-L,n-l) 
G ! * g (m,n-2) 

m.n-i 


(C.7) 


|g (ro-L,n-2) 


lg(m-L,n-L-l) 


• » 
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The  elements  bjl)  are  given  by  the  algorithm 
for  J * 1 to  L + 1 
For  I * 1 to  L 
Let  K » I + (J-l) (L+l) 

If  J ■ 1,  b!jp-  _vi+i,j 

If  j - !•  1 

If  J>1,  b».  -ibIrJ+1. -IvI+1 
if  1,  ^ i 

Otherwise  b^  = 0 
J K 

The  elements  b are  given  by  the  algorithm 

For  J * 1 to  L 
For  I * 1 to  L + 1 
Let  K * I + (J-l) (L+l) 

If  I - 1. 

(2) 

If  1 * l*bK+L,K-  1*0 

(2)  1 1 „ 

If  I>1,  b^  -Ibj-i^j  * “ J VI,J+1 

(2)  1 
If  I >1#  bK+L>K  * J 

Otherwise  * 0 

Define  the  (L+l)?  by  (L+l)  ^ matrix  A with  elements  . Here  we  depart 
frar  the  standard  notation  to  avoid  confusion  between  the  elements  of  the 
matrix  A and  the  coefficients  of  the  filter  specified  by  (C.l) . Then 

A - [aJR]  (C.ll) 


- V 


I,J+1 
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The  elements  aJKare  given  by  the  algorithm: 
For  I * 1 to  L + 1 
For  J = 1 to  L + 1 
Let  K = I + (J-l) (L+l) 

al,K  = ai-l,J-l 
Otherwise  a,  „ * 0 

J a 


Example  C.l 

Consider  the  filter  specified  by  (C.l)  or  (C.l)  where  L is  equal  to  2. 

In  that  case,  G , G . #G  and  F are  9x1  vectors  and  B,  , B-,  , and  A are 

m,n  Tn-i,n  rti,n-l  m,n  1 * 

9x9  matrices.  The  vectors  G , G , , G - and  F_  „ are  given  by 

m,n  m-l,n  m,n-I  m,n 


n 


g(m,n) 

g(m-l,n) 

g(m-2,n) 

g(m,n-l) 

g(m-l,n-l) 

g(m-2,n-l) 

g(m,n-2) 

g(m-l,n-2) 

g (m-2,n-2) 


(C.12) 
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I 

I 

l 

1 


m-l,n 


g(m-l,n) 

g(m-2,n) 

g(ro-3,n) 

g(m-l,n-l) 

g(m-2,n-l) 

g(m-3,n-l) 

g(m-l,n-2) 

g(m-2,n-2) 

g(m-3,n-2) 

• 


Sn,n-1 


g(m,n-l) 

g(m-l,n-l) 

g(m-2,n-l) 

g (m,n-2) 

g(m-l,n-2) 

g(m-2,n-2) 

g(m,n-3) 

g(m-l,n-3) 

g (m-2,n-3) 


(C.13) 


(C.14) 
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Note  that  if  we  write  the  equations  corresponding  to  the  rows  of  (C.19) 
we  obtain  for  the  first  row 


g(m,n) 


L L 

-l  T. 

J-0  K-0 


L L 

bJKg (m-J, n-K)  ♦E  I aJRf (m-J, n-K)  (C.20) 

J-0  K-0 
J+K>0 


For  the  subsequent  rows,  we  obtain 

g (m-J, n-K)  - g (m-J, n-K)  (C.21) 

or  the  outputs  are  equated  to  themselves.  It  follows  directly  that 
(C.19)  is  equivalent  to  (C.2). 

C.3  The  S Matrix 

We  now  give  the  algorithm  for  the  matrix  which  is  used  to  reorder  the 
rows  and  columns  of  the  matrices  B ^ and  B 2 for  the  proof  of  Theorem 
4 js  described  in  Appendix  A.  Define  the  (L+l)2  by  (L+l)2  matrix  S such  that 

S *[s  Jk]  (C.22) 
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The  elements  SJK  ace  given  by  the  algorithm 
For  I » 1 to  L + 1 
For  M ■ 1 to  L + 1 
Let  J - M + (1-1)  (L+l) 

Let  K » I + (M— 1)  (L+l) 

SJK“  1*0 
Otherwise  s * 0 


Example  C.2 

Consider  the  filter  specified  by  (C.l)  or  (C.2)  where  L is  equal  to  2. 
In  that  case,  we  have 


1 . 


1 1 


S 


1 0 0 0 0 
0 0 0 1 0 
0 0 0 0 0 
0 10  0 0 
0 0 0 0 1 
0 0 0 0 0 
0 0 10  0 
0 0 0 0 0 
0 0 0 0 0 


0 0 
0 0 
0 1 
0 0 
0 0 
0 0 
0 0 
1 0 
0 0 


0 0 
0 0 
0 0 
0 0 
0 0 
1 0 
0 0 
0 0 
0 1 
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A new  approach  to  the  stability  problem  for  the  two-dimensional  digital 
recursive  filter  is  presented.  The  bivariate  difference  equation 
representation  of  the  two-dimensional  recursive  filter  is  converted  to  a 
multi-input  multi-output  (MIMO)  system  similar  to  the  state  space  representation 
of  the  one  dimensional  digital  recursive  filter.  In  this  paper,  a pseudo-state  space 
representation  is  used  and  three  two-dimensional  polynomial  matrices  are 
obtained.  A general  theorem  for  stability  of  two-dimensional  digital  recursive 
filters  is  derived  and  a very  useful  theorem  is  presented  which  expresses 
sufficient  requirements  for  instability  in  terms  of  the  spectral  radii  of 
these  matrices. 


I . Introduction 
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A two-dimensional  digital  recursive  filter  can  be  characterized  by  the 

bivariate  difference  equation 

L L L L 

g(m,n)  - E 2 a f(m-J,n-K)  - E z b a(m-J,n-K)  (1) 

J-0  K-0  JK  j=0  k-0  JK 

where  the  coefficients  a and  b are  constants  [1]  and  some  of  these  constants 

JK  JK 

may  be  zero.  There  are  two  major  problems  to  consider  in  the  design  of  recursive 
filters  for  two-dimensional  signal  processing:  synthesis  and  stability.  The  synthesis 
problem  consists  of  determining  the  filter  coefficients  so  that  the  required 
frequency  response  is  realized.  If  the  resulting  filter  is  to  be  useful,  it 
must  be  bounded  input-bounded  output  (BIBO)  stable.  In  this  paper  the  stability 
problem  is  considered  and  a new  approach  to  stability  analysis  for  the  two- 
dimensional  digital  recursive  filter  is  presented. 

For  the  one-dimensional  case,  there  are  essentially  two  methods  of 
determining  necassary  and  sufficient  conditions  for  stability  of  digital  filters: 
examining  regions  of  analyticity  for  the  characteristic  polynomial  and  by 
direct  evaluation  of  the  characteristics  of  the  impulse  response  (2,3,4]. 

In  particvlar,  if  the  system  corresponding  to  the  digital  filter  is  represented 
by  a state  space  equation,  then  one  can  determine  stability  from  the  coefficient 
matrices  in  the  state  space  equation[4].  For  the  two-dimensional  case, 
generalizations  of  the  first  method  involves  examining  regions  of 
analyticity  for  bivariate  polynomials  which  is  computationally  feasible  only 
for  very  single  filters [5].  This  paper  attempts  to  generalize  the  second 
method  for  the  two-dimensional  case,  i.e.  to  establish  stability  by  computing 
the  spectral  radii  of  coefficient  matrices  with  real  coefficients. 
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II.  Pseudo  State  Space  Representation 


Fornasini  and  Marchesini  [6]  have  defined  a state  space  representation 
of  the  two-dimensional  digital  recursive  filter.  In  this  paper,  we  use  a 
particular  case  of  the  Fornasini-Marchesini  model  where  one  of  the  coefficient 
matrices  is  the  null  matrix.  Thus,  we  obtain  the  pseudo  state  space  representation 


B,  G,+  B-G  . + AF 
l m-l,n  i m,n-l  m,n 


(2) 


m,n 

g(m,n)  - DG^ 

n is  a coluun  vector  such  that  its  elements  are  the  outputs,  g(nKJ,n-K) 
where  0£J<p  and  BiKiL.  Note  that  Gm>n  contains  all  of  the  outputs  that  are 
represented  in  (1)  including  g(m,n).  Similarly,  Fm>nis  a colunn  vector 
such  that  its  elements  are  the  inputs,  f(m-J,n-K)  where  0<J<L  and  0<K_fI*. 

We  can  then  define  matrices  , B2  and  A [7]  such  that  (1)  and  (2)  are 
equivalent.  The  matrices  B^  , B2  and  A are  all  of  order  (L+l)2  by  (L+l)^. 
The  vector  D is  a row  vector  with  L+l  elements. 


The  ordering  of  the  outputs  in  G and  of  the  inputs  in  F_  is 

m,n  m,n 


not  unique.  However,  the  ordering  does  affect  the  relative  position  of  the 
elements  of  the  corresponding  coefficient  matrices.  Also  note  that  there  are 
identical  elements  in  >nand  G^n-t  Where  this  occurs,  the  corresponding 


elements  of  and  B2  can  be  divided  such  that  the  magnitude  of  each  is  no 
larger  than  that  of  the  corresponding  b jK  or  one  as  appropriate.  It 
is  convenient  to  consistently  divide  equally  and  choose  a particular  ordering 


a I 
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III.  Stability  Analysis 

The  stability  analysis  herein  will  be  confined  to  the  linear 

shift  invariant  (LSI)  two-dimensional  discrete  system.  Such  a system  is  BIBO 

stable  if  and  only  if  the  discrete  inpulse  response  of  the  system, 

h(m,n),  is  absolutely  simnable,  i.e.,  “ jh(m,n)|«®  [1]. 

m,n“0 

Let  us  define  the  particular  vector  Hj  R as  that  input  vector  which 
represents  a single  unit  sample  at  the  (J,K)  position  of  the  two-dimensional 
data  array  and  all  other  inputs  are  zero.  Let  us  further  define  the  initial 
condition  vector,  Gj.i^^and  Gj  » as  null  vectors.  Then  for  m*J  and  rpK , 
(2)  reduces  to 


GJ,K 

h(J,K)  * DGj,k  (3) 

Define  the  term  CfB^.B^)  as  the  sum  of  all  unique  products  involving  B 
as  a factor  J times  and  B2  as  a factor  K times.  It  is  helpful  to  note  that 
if  Biand  B2ccmmute,  then  C(Bi,b£)  » bJbJ'  (J+K)  !^B^/(J!K1) . In 

general,  the  matrices  do  not  commute.  Therefore,  we  give  as  an  example 
C (B  ^,B  2*  * B^B  ®2®1* 

Lenina  1:  Given  the  discrete  LSI  system  represented  by  (2) , the  contri- 
bution to  the  output  vector,  G m#n , by  a single  input  vector,  , which 
corresponds  to  a unit  impulse  at  the  (J,K)  position  where  J<M  and  K<N  is  given 
by  Gm,if  c(B^m  J»32.n  K)AHJ,K* 

The  proof  of  Lenina  1 is  given  in  the  Appendix.  Lenina  1 provides  a 
convenient  means  of  finding  the  output  of  the  two-d iemns ional  digital  recursive 
filter  for  all  values  of  m and  n when  the  filter  is  excited  by  a single 
input  at  any  point  in  the  array.  Since  the  filter  is  linear  and  shift 
invariant,  we  can  use  the  principle  of  superposition  to  find  the  output  for 


any  particular  sequence  of  inputs. 

Thus,  the  unit  impulse  response  of  the  filter  is  given  by 
C(Br,B2n)AH0f0 

h(®,n)  - DGm#n-  DC(B1m,B2n)AI^j^  0 (4) 

Lenina  2:  Given  the  discrete  LSI  system  represented  by  (2)  for  which 
the  corresponding  transfer  function  has  mutually  prime  numerator  and 
denominator  polynomials.  If  the  contribution  to  the  output  vector  G by 

m t n 

a bounded  sequence  of  input  vectors  K where  0SJ<M  and  0<K<N  can  be 

expressed  by  G - Q1^-  _or  G = „ , then  the  system  is  unstable  if 

nvn  j,k  m,n  J , K 

p(Q) , the  spectral  radius  of  Q,  is  greater  than  one.  The  proof  of  Laima  2 
is  given  in  the  Appendix. 

Theorem  1:  The  discrete  LSI  system  represented  by  (2)  is  stable  if  and 
only  if  for  at  least  one  matrix  norm 


Theorem  1 follows  directly  from  (4)  and  the  requirement  that  the 
discrete  iiqpulse  response  be  absolutely  summable.  Since  h(m,n)  is  a scaler, 

its  matrix  norm  is  equivalent  to  its  absolute  value  and  the  proof  of  theorem  1 is 
obvious. 

Theorem  2:  The  discrete  LSI  system  represented  by  (2)  and  for  which  the 
numerator  and  denominator  polynomials  of  the  corresponding  transfer 
function  are  mutually  prime  is  unstable  if  any  one  of  the  spectral  radii 
p(B^) , pOJj)  or  p (B  J+B2)  ^ greater  than  or  equal  to  one.  The  proof  of  Theorem 
2 is  given  in  the  Appendix. 

In  the  practical  application  of  two-dimensional  digital  recursive  filters, 


any  filter  with  p(Bj),  p(B2)  or  p(Bi+B2)  equal  to  one  can  be  considered  to  be 
unstable  and  should  be  avoided [8].  Goodman  [5]  has  shown  by  clever 

examples  that  two-dimensional  filters  with  nonessential  singularities  of  the 
second  kind  on  the  unit  bidisc  may  be  stable.  Such  a filter  may  have  p(Bi), 
p (B2)  or  p(Bj+B2^  equal  to  one.  However,  roundoff  errors  and  coefficient 
truncation  would  prevent  satisfactory  performance  by  such  a filter  for  most 
applications. 

Several  other  theorems  relating  to  sufficient  conditions  for  stability 
have  been  found  [7],  However,  it  has  been  shown  that  these  constraints  are 
too  restrictive  for  general  use.  That  is,  useful  stable  filters  can  be  found 
which  do  not  satisfy  the  corresponding  sufficient  conditions  for  stability. 

Computer  algorithms  are  readily  available  to  find  the  spectral  radius 
of  a matrix  with  real  coefficients.  Thus,  Theorem  2 presents  a convenient  and 
easily  implemented  technique  to  assess  the  stability  of  two-dimensional 
digital  recursive  filters. 
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APPENDIX 

In  this  appendix,  the  proofs  for  Leninas  1 and  2 and  Theorem  2 are  given. 
Since  all  vector  norms  are  equivalent,  any  convenient  norm  may  be  used  for 
either  input  or  output. 

Al.  Proof  of  Lenina  1 

We  proceed  with  a proof  by  induction.  If  we  use  (2)  and  (3)  to  obtain 
cJ+l,K'GJJt.fnd  GJ+1,-K+1  £or  input  vector  HJ,X  and  ail  initial 
condition  vectors  are  null  vectors,  we  obtain 

GJ+1,K  + b1GJ,K  = b1ahJ,K 

GJ,K+1  * B2GJ,K  * B2AHJ,K 

GJ+1,K+1  " B1GH,K+1  + B2GK+1xJ  * (B1B2  + B2Bi)AHj,K 


If  we  use  Laima  1,  we  obtain 

GJ+1,K  - C<B1'B2>ABJ,K  - B1Mj,k 
GJ,K+1  " C(Bf'Bi)AHJ,K  " B2AB j , K 
GJ+1,K+1  - C(B2'Bi)flHJiK  ■ (B1B2  + B2Bl>AHJ<K 


(A2) 


Thus  for  any  arbitrary  m and  n such  that  msj  and  n>K,  we  can  use  (2)  to 
write 

Gra+l,n  ' BlGm.n  + (S3) 

Then  using  (4)  to  find  expressions  for  Gmfnand  ^+1^-1'  we  ^ave 

Gm+l,n  “ IB1C(bJ"J,b5“J)  + R (A4) 


..  r*  ' V'. 


I ’ * * - 
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Consider  the  term,  C(Bi,E^).  All  of  the  products  in  the  term  either  have 
as  the  first  factor  or  B2as  the  first  factor.  If  Bi  is  the  first  factor 
we  must  postmultiply  by  the  sum  of  all  possible  products  such  that 
the  power  of  B^  is  decreased  by  one.  If  B2  occurs  as  the  first  factor,  we  must 
postmultiply  by  the  sum  all  possible  products  such  that  the  power  of  B 2 
if  decreased  by  one.  We  conclude  that 

C(bJ,bK)  = BjCteJ”1^)  + B2C(B1,bJ“1),  (A5) 

for  all  J and  K such  that  both  J and  K are  greater  than  or  equal  to  one.  It 
follows  directly  that 

Gm+l,„  " C(BT1-J.BrK)AHJ,K  • (A6) 

Similarly  from  (2)  we  write 

G ..  =»  B.G  . + B,G  _ . (A7) 

m,n+l  1 m-l,n  2 m,n 

Using  (4)  to  find  expressions  for  Gm_]_ t n+i  and  Gmfn»we  have 

Gm,n+1  * [B1C(Bi"J“1,b5+1“K)  + B2C(b5»_J,b5"K)]AHJ/K.  (A8) 


It  follows  that 


G . » C(B™“j,b9+1“K)AH 
m,n+l  — 1 i j#k 

Finally,  from  (2)  we  obtain 

Gm+l,n+l  “ BlGm,n+l  + B2Gm+i,n 
Using  Lemma  1 to  express  Gm  n+1  and  Gjn+1  n we  obtain 
®m+l,n+l  * + B2  C(Bf"1-J,B5-K)]AHJ>K 


It  follows  from  (A5)  and  (All)  that 

®m+X , n+1  ' c<Bf  1"J-E2+1‘K>»hj,K 


(A10) 


(All) 


(A12) 
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and  Lemma  1 holds. 


A2.  Proof  of  Lenina  2 

In  the  proof  of  Lemma  2,  we  shall  show  that  if  the  response  to  a 
particular  sequence  of  input  vectors  can  be  represented  as  given  in 
Lemma  2,  then  the  system  is  unstable  if  pCQpl  [91. 

Define  the  eigenvalue  corresponding  to  the  spectral  radius  of  Q as 
and  the  corresponding  eigenvector  as  1^.  Then  if  the  system  transfer 
function  has  mutually  prime  nunerator  and  denominator  polynomials  we  can 
select  an  input  vector  such  that 


AF  » tPn  +R_  _for  all  J and  K 
J,K  Q J,K 


(A13) 


where  e is  an  arbitrary  nonzero  finite  constant  and  RJK  is  not 
in  the  direction  of  Pq.  We  then  have 

Vn  ’ '<^0  + 0" RJ,K  <u4> 

'Rjen  since  Xq  is  the  eigenvalue  corresponding  to  the  spectral  radius, 

the  norm  of  n is  dominated  by  the  term  £QPq  in  the  limit  as  m approaches 
infinity. 

Thus 


S 


Lim 

m-**» 


Note  that  S is  infinite  if Xq  is  greater  than  one  and  Lenina  2 holds. 

A3.  Proof  of  Theorem  2 

For  this  proof,  we  show  that  we  can  find  a particular  sequence  of  inputs 
that  give  unbounded  output  if  either  of  the  spectral  radii  specified  in 
Theorem  2 is  greater  than  one. 


From  Leona  1 and  2 the  output  from  a single  arbitrary  bounded  input  at 
the  (J,K)  position  can  be  given  by 

GM,N"  (A16) 

g(M,N)  - DGM/N 

vAiere  F(J,K)  is  the  scalar  input  at  the  (J,K)  position.  If  we  let  K*N  and 
J=0  in  (A16) , we  have 

gM,N*  F(0,N)C(Bi,B§AHj  f f(0,N)^RBj#K  (A17) 

If  we  apply  Lemma  2,  we  see  that  the  system  is  unstable  p <Bi)>l.  If  we  let 
J * M and  K * 0 in  (A16) , we  have 


GM,N  -*(M.0)C(b“,B^)AHj  k - f(M,0)B>J  K 


(A18) 


If  we  apply  Lenina  2,  we  see  that  the  system  is  unstable  if  p (B2)  >1. 

If  we  use  a particular  sequence  of  inputs  f(J,M-J)  for  where  all 

f(J,M-J)  are  bounded  and  equal.  Using  the  principle  of  superposition 
and  (A16)  we  have 


G»'1<  ' L 


(A19) 


Since  all  inputs  are  equal,  we  can  write 

gm,n  ■ ‘to  c(bi’J'b?,1aho,K 


- f (0,M)  (B.  + B0)MAH 


M,M 


0,M 


(A20) 

(A21) 


If  we  apply  Lenina  2,  we  see  that  the  system  is  unstable  if  p (B1-*fi2)>1* 
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INTRODUCTION 


Technology  would  be  much  advantaged  by  the  availability  of 
better  methods  for  anticipating  and  preventing  the  mechanical 
failure  of  ceramics.  In  particular,  the  dependence  of  mechanical 
properties  on  the  microstructural  character  is  so  complex, 
and  ceramic  microstructures  so  inherently  variable,  that  the 
technological  use  of  ceramics  has  been  persistently  retarded. 

The  field  of  mechanical  properties  is  vast,  therefore  this 
phase  addressed  itself  to  a well  defined  and  limited  objective; 
to  acquire,  through  literature  reviews  and  laboratory  testing, 
and  understanding  of  the  principles  of  microhardness  testing  of 
ceramics.  This  understanding  was  then  used  to  devise  a modification 
of  the  microhardness  test  which  allowed  the  extraction 
of  additional  information. 

The  isolation  of  the  hardness  test  as  the  centerpiece  for  this 
task  can  be  easily  defended.  Hardness  tests  have  long  been  used  for 
classification  and  surveying  of  various  materials,  and  their 
application  to  ceramics  promises  to  be  useful,  for  example  to 
determine  fracture  toughness.  And  in  a general  way,  indentation  serves 
as  a model  system  for  studying  the  strength  degradation  in  a wide 
range  of  phenomena  such  as  machining  flaws  and  impact  abrasion. 


HARDNESS  TESTING 


The  hardness  test  has  long  suggested  itself  as  a simple  means  of 
obtaining  mechanical  properties:  the  ease  of  the  test  and  the  snail 
sample  size  needed  being  the  principal  inducements.  As  a result,  a 
number  of  technically  valuable  correlations  have  been  established, 
for  example,  the  tensile  strength-Brinell  hardness  formula  for  certain 
metals.  Unfortunately,  the  hardness  test  has  not  become  a routine  tool  for 
investigating  ceramic  materials.  The  principal  factors  in  this  was 
the  controversy  in  the  relative  roles  of  plastic  deformation  and 
cracking,  (it  was  not  until  the  70's  that  plasticity  was  accepted  as 
an  important  process) , and  the  lack  of  theoretical  underpinning  to 
sort  out  confusing  results. 

The  mechanics  of  point  indentation  has  been  slowly  developing. 

As  early  as  1881  Her  tz^  analyzed  the  gene  rad  elastic 
contact  between  two  curved  bodies,  and  in  1885  Boussinesq^2)  solved 
the  stress  field  for  the  case  of  an  infinitely  sharp  indentor  on  a 
flat  surface.  However,  as  the  stress  field  produced  in  non-homogeneous, 
anisotropic  crystalline  materials  is  complex,  and  there  are  no 
solutions  for  stress  fields  in  terms  of  general  geometry  for  an 
elastic  indentor,  it  seems  unlikely  that  any  solution  will  emerge. 

Anyone  who  has  made  hardness  tests  on  ceramics  will  have  been  impressed 
with  the  diversity  of  results  which  arise  from  minor  variation*? 
in  test  conditions. 

Nevertheless,  a number  of  models  have  been  devised  to  understand 
the  indentation  of  actual  materials,  containing  drastic  simplications 
and  based  on  systematic  studies  of  ceramic  materials^3*4)  . Most  of 
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these  models  have  in  cannon  the  inclusion  of  elastic-plastic  processes, 
and  the  assunption  of  spherical  symmetry;  in  addition,  the  initiation 
and  propagation  of  cracks  are  usually  treated  as  separate  events.  There 
are  also  two  extreme  categories  of  contact  situations;  blunt  and  sharp. 

For  blunt  indentors*5}  the  crack  nucleates  from  pre-existing  flaws  and 
develops  into  a Hertzian  cone,  while  sharp  indentors^)  nucleate  cracks 
from  the  plastic  zone  at  the  contact,  which  then  develops  into 
half-penny  cracks. 

Recent  theoretical  and  experimental  results  have  attempted  to 
establish  the  actual  macroscopic  events  in  indentation  fracture.  The 
post-mortem  examination  of  fracture  surfaces  and  indentation  impressions 
have  been  important  in  directing  theoretical  attempts.  Perhaps  the 
most  dramatic  development  is  the  apparent  correlation  between 
fracture  mirror  patterns  and  fracture  stress  for  glass  and  ceramics^7-12). 
Unfortunately,  after  much  intense  study  and  confirmation  of  the  general 
empirical  relationships,  it  appears  that  these  relationships  may  have  no 
useful  fundamental  implications ^ 13-15) . 

Not  withstanding  such  analytical  limitations,  a number  of  studies 
have  established  certain  common  features.  The  greatest  attention  has  been 
directed  to  the  blunt  indentor;  the  crack  systems,  their  nucleation, 
propagation,  and  geometry  have  been  studied . Generally,  the  blunt 
indentor  probes  the  cleavage  tendencies.  The  effect  of  loads^22*  and 
indentor  angle  (sharpness) (3»6> have  been  studied.  The  theoretical  under- 
standing of  the  crack  growth  in  these  tests  has  been  primarily  based  on  the 
Griffith  energy  criterion. 

The  sharp  indentor  is  probably  more  pertinent  to  actual  contact 

wm  LJ 

E—  - 


,9 


situations,  such  as  grinding  and  hardness  testing.  The  principal  difference  from 
the  blunt  indentor  is  the  plastic  flow  preceeding  the  formation  of  the 
crack  and  the  complex  stress  field.  Once  the  crack  system  develops,  the 
influence  of  indentor  geometry  becomes  less  important.  Specific  discussion 
of  the  sharp  indentation  studies  will  be  in  terms  of  the  hardness  test,  which 
will  be  considered  a semi-sharp  indentor. 

In  stannary,  while  there  have  been  determined  various  functional 
relationships  for  indentation  processes,  the  theoretical  description 
remains  incomplete. 

HARDNESS  TESTING  OF  CERAMICS 

In  spite  of  the  formidable  theoretical  and  experimental  difficulties 
a number  of  attempts  have  been  made  to  use  the  hardness  test  on  ceramics. 

Most  of  these  have  involved  the  port-mortem  examination  of  crack  patterns 
and  dimensions,  and  the  indentor  has  most  frequently  been  the  Vickers 
diamond  pyramid.  This  indentor  is  usually  considered  a semi-sharp  indentor 
and  its  indentation  stress  field  is  favorable  for  plastic  flow.  Specific 
examples  sure  the  measurement  of  stresses  in  tempered  glass  surfaces^ 23 ^ . 
fracture  toughness  determination ^24) , and  compressive  strengths^23* . 

The  results  of  some  of  these  will  be  referred  to  in  the  discussion  of 
experimental  results. 

EXPERIMENTAL  CONJECTURES 

The  experimental  idea  developed  here  is  an  extension  of  the  attempts 
to  use  the  features  of  the  hardness  impression  such  as  the  extent  and 
pattern  of  cracking.  The  surface  crack  pattern  as  a clue  to  fracture 
toughness  was  originated  by  Palmquist  in  1957 (27) and  since  has 
been  considerably  refined  in  theory  and  procedure ^ 28-29 ) . The  shortcoming 


of  these  methods  is  that  measurement  of  the  crack  features  is  difficult 
and  unreliable,  usually  requiring  SEM  and  etching  techniques* 30-3  2 > . 

A sinplier  and  less  ambiguous  method  is  desirable. 

Although  there  is  no  solution  of  the  stress  field  for  hardness  indentors, 
the  deformation- fracture  for  well  developed  cracks  in  brittle  materials 
is  known,  and  can  be  used  to  orient  our  thoughts.  Fig.  1 shows  a Vickers 
Diamond  Pyramid  Indentor  (DPI)  and  the  primary  crack  pattern.  The  sharp 
edges  of  the  indentor  tend  to  initiate  aind  develop  half-penny  cracks 
along  the  diagonals.  Therefore,  to  prove  the  extent  of  these  cracks  it 
seems  reasonable  to  search  for  some  interaction  of  the  cracks  with  stress 
singularities  such  as  a surface  or  another  indentation.  For  example,  as 
shown  in  Fig.  2,  for  the  orientation  of  indentor  shown,  the  cracks  would 
extend  to  the  surface,  and  the  crack  length  would  be  unambiguously 
revealed  when  breakthrough  at  the  edge  occurred.  Another  possibility 
is  shown  in  Fig.  3,  where  the  distance  between  successive  impressions  is 
varied  until  chipping  occurs. 


Figure  1.  Diamond  Pyramid  Indentation  and  Resulting 

Crack  Pattern 
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EXPERIMENTAL  MATERIAL 

The  materials  used  were  remnants  from  an  evaluation  program  on 
ceramic  vane  materials  conducted  by  I IT  Research  Institute  and  the  Air 
Force  Materials  Laboratory  and  were  extensively  characterized  by  theirf33"35>  . 
Hie  materials  are  briefly  described  below,  and  some  properties  summarized 
in  Table  1.  For  testing  by  North  Carolina  A & T State  University,  the 
surfaces  were  polished  with  diamond  paste,  final  polishing  being  with 
1 -y  diamond. 

NC-350:  Reaction  bonded  Sij  . The  microstructure  has  25%  porosity, 
with  the  porosity  and  silicon  nitride  phase  uniformly  distributed.  The 
porosity  was  one-half  to  one-third  open.  The  phases  present  were 
« -si3N4  (major)  and  6 -si3N4 (minor) . 

NC-435:  Siliconized  SiC.  A two-phase  material  with  about  20% 

Silicon.  The  SiC  phase  is  the  « -form,  and  low  porosity  which  is 
mostly  closed. 

NC-132:  Hot  pressed  Si^4.  The  phases  present  were  « -Si3N4,  the  major 
phase,  and  Si2N40,  the  minor  phase.  There  were  also  traces  of  WC. 

The  microstructure  showed  fine,  elongated  grains  with  virtually  no 
open  porosity. 

EXPERIMENTAL  APPARATUS 

The  tests  urere  conducted  on  an  Kentron  microhardness  tester,  using 
a 136*  Diamond  Pyramid  Indentor.  This  indentor  is  cut  in  the  shape  of  a 
square-based  pyramid  having  an  apex  angle  of  136*.  The  loads  were  applied  for 
15  seconds,  and  care  was  taken  to  reduce  machine  vibration.  The  tests  were 
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j.  Table  1 


Material 

Batch 

Eracture  Stres 

& Modulus 

tensity 

NC-350 

1 

30,860  psi 

25.5x10  psi 

2.52  gm/cc 

2 

23,530 

23.5 

2.40 

3 

33,960 

27.7 

2.54 

NC-435 

1 

50,460 

53.9 

2.936 

3 

65,980 

49.5 

2.997 

4 

55,150 

48.8 

2.962 

NC-132 

1 

90,910 

48.4 

3.177 

2 

115,210 

45.7 

3.186 
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RESULTS 

Based  on  the  previous  speculations,  a nunber  of  different  indentation 
procedures  were  tried.  The  two-indentor  method  did  yield  seme  success; 
however,  the  experimental  procedure  proved  troublesome.  At  least  on  the 
machine  used,  it  proved  difficult  to  line  up  the  impression  accurately, 
and  so  not  enough  usable  data  resulted.  The  test  using  an  edge,  shown  in 
Fig.  3,  gave  better  results.  In  this  test,  the  distance  from  an  edge  which 
caused  breaking  of  the  edge  at  a particular  load,  was  proved.  The  results 
are  given  in  Tables  2 and  3.  No  results  are  included  for  NC-132,  as  this 
material  frequently  deformed  by  flaking  rather  than  forming  an  impression 
and  cracking.  A few  usable  data  points  indicated  that  the  results  for 
this  material  followed  tht  trends  set  by  the  other  materials.  This 
material  is  included  because  it  apparently  represents  the  upper  strength 
limits  for  materials  whish  can  be  used  for  the  distance  from  the  edge  test, 
although  anisotropy  might  be  playing  a role. 

DISCUSSION 

The  data  of  Table  2 are  plotted  in  Pig.  4 is  the  fracture  stress 
as  determined  by  four  point  bending  from  Table  1.  While  there  is  no  single 
straight  line  which  appears  to  fit  this  data  well,  it  is  possible  that  a 
line  can  be  fitted  through  each  of  the  sets  of  points  for  the  two 
materials.  Certainly  the  hardness  test  appeared  to  rank  the  materials 
correctly,  both  as  to  batch  and  type;  this  certainly  indicates  acme  accuracy 
and  shows  that  surface  finish,  machine  operation  and  microstructural 
variability  are  not  scattering  the  data  beyond  use. 

The  data  from  Table  3 are  shown  plotted  on  log-log  paper  in  Fig.  5. 
Comparison  of  this  data  with  that  in  Table  1 shows  that  the  materials  are 


Table  2 


Material 

Batch 

■Hardness* 

NC-350 

1 

859 

2 

662 

3 

1087 

NC-435 

1 

1086 

3 

1780 

4 

1646 

600  gram  load. 


♦Average  of  10  readings, 
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ranked  by  modulus  rather  than  by  strength.  In  trying  to  understand 
this  and  other  features  about  this  data,  we  need  some  theoretical 
underpinning.  The  first  obvious  correlation  to  look  for  is  some 
analogy  to  the  mirror  boundary  relationship (9>10,14) 


where  V.  is  the  fracture  stress,  r the  fracture  mirror  radius, 
r in 

and  A the  mirror  constant.  Although,  as  mentioned  previously,  the 
theoretical  basis  for  this  relationship  has  been  denied,  it  has 
nevertheless  been  confirmed  by  many  studies,  and  so  there  is  some 
expection  that  the  present  data  would  fit  this  equation.  If  the 
distance  from  the  edge  is  taken  as  an  indication  of  the  mirror  size, 
then  the  functional  dependence  in  the  above  equation  must  be  reversed. 
That  is,  the  mirror  boundary  represents  the  surface  created  by  the 
fracture  stress,  so  that  a strong  material  vrould  create  a large 
fracture,  while  in  the  present  test,  as  the  load  is  constant  (or 
normalized) , the  stronger  material  will  break  closer  to  the  edge  and 
thus  create  a smaller  area. 

Even  though  the  functional  dependence  is  correct,  the  data 
slope  is  not  the  required  1/2,  and  no  amount  of  manipulation, 
for  converting  the  distances  to  areas,  can  make  it  conform  to  the 
mirror  boundary  equation. 

A second  approach  is  to  consider  the  distance  from  the  edge  as 
proportional  to  the  characteristic  crack  dimension.  This  would  mean  that 
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the  crack  extension  stage  would  predominate,  and  that  the  force  would  be 

proportional  to  the  3/2  power  of  the  crack  dimension  (or  distance  from 

the  edge),  as  expected  for  a half-penny  crack (34,35) # The  data  does 

appear  to  cluster  around  a 3/2  slope,  although  it  is  impossible  to  know 

if  the  scatter  in  slopes  is  due  to  test  or  material  variability,  or  is  a 

natural  material  function.  If  the  distance  from  the  edge  is  truly  an 

indicator  of  the  force  necessary  to  extend  a crack,  then  this  method  may 

eventually  make  it  possible  to  determine  such  factors  as  the  fracture 

surface  energy.  IITRI  determined  K by  the  double  torsion  method  for 

C 

the  materials  tested  here.  See  Table  4.  Those  values  correlate  to  the 
ranking  determined  here,  if  the  small  differences  between  batches  3 and  1 
are  ignored. 

If  we  then  review  all  of  the  data,  the  strengths  of  these  mater ials 
correlate  with  the  densities,  which  correlate  with  the  hardness, 
and  the  modulus  correlates  with  the  K which  correlates  with  the  distance 
from  the  edge.  Also,  a plot  of  the  distance  from  the  edge  vs  load  gives 
roughly  the  slope  3/2  for  all  batches,  the  relationship  between  the  force 
and  extension  of  a half-penny  crack.  Intriguing  as  these  results  are, 
it  vould  be  overopt imistic  at  this  time  to  suggest  that  anything  other 
than  a method  for  ranking  materials  or  batches,  by  modulus  or  K , has 
been  developed.  An  extremely  limited  range  of  materials  was  tested,  and 
as  the  difficulties  with  the  NC-132  material  showed,  some  materials  may  be 
difficult  to  accomodate  to  this  test.  It  is  hoped  that  these  results 
will  encourage  further  exploitation  of  the  hardness  test. 
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INTRODUCTICN 

The  second  project  of  this  investigation  was  entitled  "Strength 
characterization  of  brittle  materials  by  means  of  s imple— model  destructive 
testing  and  by  indentor- initiated  crack  growth". 

A simple  destructive  test  is  appealing  due  to:  (1)  its  decisiveness  in 
quantifying  failure;  (2)  its  ability,  if  successful,  in  transferrability  of 
results  to  complex  designs  and  shapes;  (3)  its  incorporation  of  size  effects; 
and  (4)  in  that  it  requires  inexpensive  testing  equipment. 

Literature  on  ceramics  largely  resides  in  three  camps,  fracture  mechanics, 
continuun  mechanics,  and  statistical  failure.  Fracture  mechanics  predictions 
of  failure  couple  with  elasticity  theory  where  strain  energy  release  rate  is 
equated  with  a critical  flaw  length  or  size  [1-2]  frequently  on  the  basis  of 
standard  'crack'  models  such  as  interface  delamination  in  the  case 
composite  material  problems  or  as  voids  with  sharp  edges  as  in  the  case  of 
the  penny- shaped  crack  [2].  The  difficulty  in  such  approaches  involving 
elasticity  theory  is  the  intractibility  of  the  field  equations  even  for  linear 
elastostatic  problems  of  simple  geometry  [3,4]. 

Finite  element  and  difference  methods  offer  some  advances  in  the  numerical 
techniques  [5,6],  in  particular  when  coupled  with  cine ill ary  methods  that 
predict  the  stress  intensity  factor  [7,8].  These  methods  all  fail,  however, 
to  give  an  account  of  the  customarily  larger  experimental  scatter  in  apparent  failure 
strengths  of  nunerous  grainy  brittle  materials. 

The  classical  work  on  statistical  failure  by  Weibull  [9-11]  has  been  applied 
to  nunerous  stocastic  processes  capacitor  discharge,  fatigue  time  life,  and 
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others  largely  due  to  the  ease  of  use,  simplicity  of  assumptions  and  non- 
stocastic  processes  as  ceramic  fracture  in  relation  to  a size  effect. 

The  major  problems  posed  by  this  theory  are  the  mathematical  determination 
of  the  Weibull  parameters,  in  particular  for  the  three  parameter  distribution, 
and  in  the  application  of  the  results  in  a meaningful  manner  to  the  elasticity 
solution  in  order  that  a failure  probability  for  the  body  can  be  found. 

In  order  to  develop  a technique  that  would  utilize  Weibull  Theory 
and  an  elasticity  solution  to  generate  fracture  probabilities  for  a body  of  brittle 
material,  it  was  necessary  to  investigate  what  techniques  were  proposed  and  used  to 
generate  the  Weibull  parameters  to  begin  with. 

Sample  earlier  work  was  soon  found  based  on  moment  generating  methods 
and  on  rank-order  theory  [12,13]  to  differ  with  results  obtained 
by  several  iterative  methods  whereby  the  sume  of  squares  of  deviations  were 
minimized  (residuals)  as  is  the  case  of  standard  least-squares  techniques. 

Moment  generation  can  involve  large  truncation  error  and  rank-order  assumptions 
should  be  significantly  less  precise  than  comparison  with  cumulative  totals 
of  experimental  data.  Further,  although  least-squares  analyses  tends  to 
emphasize  deviations  in  the  experimental  curve  ordinates,  incorporation  of 
appropriate  weights  can  counter  this  effect  while  preserving  the  methods 
basic  simplicity. 

Another  problem  with  earlier  techniques  is  that  an  assumption  was  generally 
required  for  multiaxial  stress  states  necessitated  by  the  uniaxial  stress 
state  under  which  data  had  been  taken.  Many  persons  assumed  statistical 
independence  in  the  actions  of  principal  direction  stresses,  i.e.,  the 
probability  of  failure  of  an  element  equal  to  the  product  of  probabilities 
of  failure  due  to  each  principal  direction  stress  [12].  Some  writers  have  proposed 


extensions  of  the  uniaxial  Weibul  Theory  based  on  fracture  mechanics  where 
distributions  of  penny-shaped  cracks  initiate  failure  [14,15,16], 

In  this  report  it  was  proposed  to  couple  the  results  of  multiaxial  tests, 
notably  biaxial  ones,  with  an  analysis  and  development  of  Weibull  parameters 
obtained  at  specific  stress  ratios.  This  would  aid  the  utilization  of  experi- 
mental data  without  the  need  to  develop  further  fracture  mechanics  assumptions. 
Further,  it  has  been  hoped  that  such  a technique  when  coupled  with  effective 
testing  of  standard  specimens  and  stress  fields  could  help  lend  a standardization 
in  the  analysis  of  experimentally  obtained  statistical  fracture. 

DEVELOPMENT  OF  SIMPLE  TEST  ANALYSES 

An  extension  to  publications  where  standard  specimen  shapes  are  used  to 
generate  Weibull  statistics  [17,18]  are  developed  below  where  the  three 
parameter  Weibull  fit  is  used. 

UNIAXIAL  CONSTANT  STRESS  FIELDS 

Hie  Weibull  distribution  for  a unaxial  constant  stress  state  is 
given  below: 

F(O)  - 1 - e x p/nc  ff°  Su'T  du 

u - Unit  Volume 
on 

K,  ou,  cr#,  m ■ weibull  Parameters  which  can  be  expressed  uniquely  as 
three  independent  ones. 

For  unaxial  constant  stress,  is  a principal  stress,  independent  of  the 
volune  so  that 


-J 'a  .. , 
vm«h)m 


Upon  using  C 
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r«J)  - 1 - axp(-c(0-<u)») 

The  simpliest  least-squares  solution  is  obtained  by  noting  that  taking 
logs  of  this  equation  twice  gives: 

Inlnf- i £n(c)  + a la(a-o) 

\i  - r(o>y  u 

This  is  a linear  relationship  allowing  a linear  least-square  analysis 
to  be  invoked  for 
y ■ A+Bx 

where  A ■ JLn(c) 

B - a 

Xx  - X(Ol)  - 
,<at,  - y,  - 

l • 1,2,...,  n - the  number  of  data  points. 

The  solution  that  minimizes  the  residual  is: 

A * ^Zwlylxl)  (Zwx)  - (Iwlzl)  (Zwlyl))/D 

B ■ ((Zw^x2)  (Zwlyl)  - (Zwlxl)  (Zw1xlyl>)  /D 

D - <Zwx)  (Zwxx^)  - <Zwxxx)2 
Where  all  iinaatioM  are  l - 1,...,  n. 

In  case  equal  weights  are  desired,  the  quantity  is  replaced  by  n, 
otherwise,  wx  » 1.  The  obvious  difficulty  with  this  analysis  is  that  A and  B are 
dependent  on  ou  . To  determine  the  best  on,  an  iterative  programing  scheme. 
Appendix  I,  was  developed  which  recomputes  the  residual  sun  of  squares 
and  chooses  the  on  minimizing  this  function.  By  a slight  modification  this 
scheme  can  be  used  to  redefine  the  weights  iteratively  to  emulate  any  curve- 
fitting power  law. 
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PURE  BENDING  STRESS  FIELD  IN  A RECTANGULAR  SECTION  BAR 

Tests  utilizing  bending  stress  distributions  are  connon  and  generally 
less  expensive  than  ceramic  tension  tests  — particularly  in  the  base  of  three 
and  four  point  loading  tests.  Following  the  two  parameter  derivation 
[17] , where  the  stress  distribution  is  variable  and  must  be  integrated 
over  the  volume  in  bending,  the  probability  of  failure  is  given  by: 


Where  V - bLh 


Figure  1.  Prismatic  beam  under  4-point  loading  and  distribution  of  extri 
fiber  a tress 


Fig.  1 depicts  this  case  with  its  test  volune.  This  distribution  follows 
the  seme  form  of  the  case  of  a constant  uniaxial  stress  field  and  may  be 
solved  interatively  by  the  technique  where  upon  taking  Inin , one  has 
Inin^j ptJJ”)) " *■“ <v/2) -An + -**n *ao * 


A - in(V/2)  - in(*fl)-m  into#) 


PURE  BENDING  FOR  A CIRCULAR  CROSS-SECTICN  ROD 


Figure  2.  Bending  for  e circular  cross-section  rod 

must  be  evaluated  as  a function  of  and  M before  the  previously  discussed 
Inin  least-square  procedure  can  be  applied.  Fig.  2 depicts  this  case  and 
its  test  volune,  VT. 


results 


PURE  TORSION  FOR  A CIRCULAR  BAR 

Wider  the  action  of  pure  torsion,  the  stress  field  t - ZE 

rz  J 

from  which  the  principal  direction  stresses  are: 

* Trz'  °2  “ • 

In  order  that  a Weibull  distribution  based  on  this  test  can  be  evaluated, 
two  parameters,  discussed  in  the  next  section  must  be  defined: 
o«-  Oj  + a\ 

and  8 m Tan-^-  (—22.)  m 135  ° 

ai 

They  are  the  intensity  of  the  biaxial  stress  field,  0 , and  the  angle  on 
the  biaxial  stress  distribution  that  the  point  on  the  fracture  surface 
makes  with  the  0!  axis,  as  depicted  below: 


Fig.  3.  Biaxial  Stress  Envelope 


82 


In  this  case,  the  stress  distribution  is  integrated  over  the  volume  as 

before  resulting  in 
P«J)  - 1 - e”Bn 

Wt«r.  B - -2li_ 

(a+l)of  J O a2  (*+2) 

This  result  has  the  same  form  as  the  case  of  the  uniaxial  tension  field  and  can 
be  correlated  with  the  three  parameter  Weibull  distribution  by  the  Appendix  I 
program. 

FOURIER  ANALYSIS  OF  BIAXIAL  DATA 

All  three  Weibull  parameters  are  functions  of  the  principal 
direction  stress  ratio  or  0 where 

9 - Tan"1 


Thus  each  Weibull  parameter  c,Oj ,m  could  be  written  by 
c - c (0) 


1 

[c 

o#-  0,  (0) 

or 

We  - . 

°b 

L“J 

m ■ (0)  ■ We(0±2nir) 

invoking  periodicity  we  also  have 
We  (6)  - We(0±2nir) 

This  same  theory  can  be  developed  for  the  solid  angle  representation  in 


T 

L 
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the  case  of  three  principal  direction  stress  components; 

that  extension  is  based  on  orthogonal  polynomials  and  representation 

completeness  of  harmonic  analysis  [3,4,19]. 

Once  the  three  Weibull  parameters  are  known  for  a variety  of  8;  , the 
equivalent  parameters  can  be  obtained  for  any  angle  by  means  of  trigonometric 
interpolation.  The  computer  program  FOURIE,  Appendix  II,  requires  the  Weibull 
set  at  a variety  of  . At  this  point,  the  highest  harmonic  obtainable  for 
the  nunber  of  data  points  is  used  for  the  interpolation  series  or  the  highest 
harmonic  is  used  which  satisfies  the  null  hypothesis  condition  [20] . The 
advantage  of  this  program  is  that  the  coefficients  to  the  Fourier  harmonics, 
once  generated,  can  be  used  for  each  element  stress  ratio  without  additional 
coefficient  computation. 

The  computer  program  vrould  naturally  be  much  simplier  with  equally-spaced 
Qx  data  points  but  such  is  not  sufficiently  general,  hence  orthogonality 
cannot  be  invoked  without  the  creation  of  special  orthogonal  sequences. 

Basically  given  n (n  odd)  angles  6t  x-i,...,n,  a set  of  functions  can  be 
generated  with  l-B."1  to  define  ya(9)  where  M L.  In  this  case: 

y (0)  « _S_  + I a*,  coa(j0)  + b“  sin(i0)  . 

* 2 ,.l'  ‘ 

Hie  residual  function  is: 


H(a?,...,aJ  - E w(0l)5VyB(0){2 
i«l  J 


This  function  is  dependent  on  the  weights  Associated  with  each  Weibull 
3-9et  we(0t)-f^  It  can  be  shown  that  evaluation  of  the  above,  minimizing  H, 
leads  to  the  partitioned  matrix  equation: 
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1 Zwl  «*(V  **••  sin(0x)  .... 

«,  IWX  008(0^  ZWl  OMS(6l)  ....  COS  (0 1 ) 810(0^ 

• • to 

• • 

* 

• • • 

«,  Iwx  sin(0x)  Ivx  costf^  sin(0x)  ....  E*x  sin2(0x)  .... 


Summations  are  \ - 1, . . . ,n. 

For  brevity  only  the  first  term  in  each  set  is  listed  here.  Modification 
of  first  row  terms  leads  to  symmetry  of  the  coefficient  matrix. 

FINITE  ELEMENT  ANALYSIS 

The  remaining  task  in  the  use  of  biaxial  Weibull  statistics  is  the  use 
of  a general  stress  analysis  method  to  find  the  stress  distribution  in  a 
body  so  that  this  theory  can  be  applied.  The  crux  of  this  process  is 
a finite  element  analysis  model  written  by  J.  Brisbane  [21]  for 
bodies  of  revolution  loaded  axisymmetrically.  This  code  has  been  modified 
to  create  files  on  the  principal  stresses  for  each  element  as  well  as 
volune  for  a particular  material  designation.  The  modified  program, 
Appendix  III,  does  generate  triaxial  stress  data  but  for  most  loading 
conditions,  one  stress  component,  in  or  near  the  r -coordinate  direction 
can  usually  be  neglected. 


Zwifi 


jZwxfxcos(0 


Pwlfxsia(0 
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As  with  all  programs  of  this  type,  input  information  is  comprehensive 
and  complex.  Only  changes  are  described  when  they  relate  to  the 
Vfeibull  analysis.  One  additional  parameter,  the  material  nunber 
N is  required,  this  is  the  nunber  of  the  material  subject  to  the  Weibull  modeling. 

In  order  to  analyse  Weibull  statistics,  two  additional  disk  files  have 
been  added,  unit  43  and  unit  47.  Unit  43  stores  the  location  of  each 
model  point  of  each  element  of  material  N.  it  also  stores  the  volume 
on  each  record.  Unit  47  records  for  each  material  N element,  the  principal 
stresses.  These  files  are  saved  after  completion  of  execution  so  that  a 
Weibull  summation  can  be  performed  in  conjunction  with  the  output 
of  Appendix  II.  The  program  of  Appendix  requires  data  from  the 
Appendix  I program  or  a modification  thereof. 

CONCLUSION 

IXiring  this  next  years  program,  this  Weibull  analysis  will  be 
performed  and  compared  to  data  derived  from  experiments  on  test  specimens. 

Each  experiment  is  designed  to  provide  Weibull  parpameters  for  specific 
principal  stress-ratios  and  will  complement  entirely  the  present 
vork.  In  order  to  compare  this  biaxial  theory  to  other  methods,  one  test 
specimen  is  being  designed  with  a variety  of  stress  ratios  present  within 
its  volume.  Its  fracture  behavior  will  be  noted  and  compared  to  predictions 
based  on  this  and  the  independence  of  principal  direction  theory  of 
failure.  Thus  vhat  has  remained  is  am  experimental  confirmation 
of  this  theory  which  should  be  provided  within  the  1978  - 1979  period. 
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C PROGRAM  WEIBUL  THREE  PARAMETERS 
C THE  PARAMETERS  ARE  SIGMAO,  M , AND  C 
C THIS  PROGRAM  COMPUTES  ALL  PARAMETERS  BY 
C LINEAR  REGRESSION  METHOD.  FOR  NOTATION  SIGMAO-SIGO 

DIMENSION  X( 20 ) , F( 20 ) , W( 20 ) , EK( 20 ) , P( 20 ) , WN( 20 ) , WM( 20 ) 
WRITE(5 , 10) 

10  FORMAT( 1H.8HG0  AHEAD) 

C GO  AHEAD  REQUIRES  YOU  TO  TYPE  NO.  OF  DATA  FILE 
READ(5 ,*)  ND 
READ(ND ,*)  N 

C ND-NUMBER  OF  DATA  FILE 
C N-NUMBER  OF  WEIBUL  POINT  PAIRS 
C W- ARBITRARY  WEIGHTS 
C X-FRACTURE  STRESS  VALUES 

READ(ND,*)  ((X(I),F(I),W(I)),I-1,N) 

DO  1 I-1,N 

1 WRITE(5 , 20)  ( I ,X(I) ,F(I) ,W(I) ) 

20  F0RMAT(1X,°X(I),F(I),W(I)  FORI-0 ,I2,3(1X,E12.6)) 

DO  2 1-1 ,N 

EK(I)-ALOG(ALOG(l ./(l.-F(I)))) 

2 CONTINUE 
NN-5 
JK-0 
INT-10 
RESO-l.E+30 
BEG— .99*X(1) 

EMD-.99*X( 1) 

98  H» ( EMD-BEG ) / FLOAT ( INT ) 

DO  95  I-O.INT 
SIK-BEG+FLOAT(I)*H 
DO  96  II-l.N 

96  P( II )-ALOG(X( II)-SIK) 

CALL  LIN(W,WN,P,EK,N,SIK,A,B,RES) 

IF(RES.GT.RESO)  GO  TO  95 

RESO-RES 

SIG-SIK 

AG-A 

BG-B 

J-I 

95  CONTINUE 

WRITE(5 ,80)  AG, BG, RESO, SIG , SIGO 
C REDIFINE  EMD  AND  BEG  AND  CONTINUE 
IF  (JJC.GT.NN)  GO  TO  97 
JK-JK+1 

IF  (J.EQ.O)  J-l 

IF  (J.EQ.INT)  J-INT-1 

EMD-BEG+FLOAT( J+l )*H 

BEG-EMD-2.*H 

RESO-RES 

SIGO-SIG 

GO  TO  98 

80  FORMAT (2X,° AG- 0 ,E12.6, 7X,#BG-° ,E12.6/ ,2X, °RESO-° ,E12 .6, 


05300 

05400 

05500 

05600 

05700 

05800 

05900 

06000 

06100 

06200 

06300 

06400 

06500 

06600 

06700 

06800 

06900 

07000 

07100 

07200 

07300 

07400 

07500 

07600 

07700 

07800 

07900 

08000 

08100 

08200 

08300 

08400 

08500 

08600 

08700 

08800 

08900 

09000 

09100 

09200 

09300 

09400 


1 5X,°SIG*° ,E12.6,5X,°SIG0*0 ,E12.6) 

97  EC-EXP(AG) 

ESIGO-SIG 

EM-BG 

WRITE (5 , 100)EC ,EM,SIG 

100  FORMAT(/ , 2X, °EC*° ,E12 .6,5X, °EM*# ,E12 .6, 5X, °SIG“° ,E12 . 6/ ) 

C CALCULATED  PROBABILITIES  ARE  LISTED  AS  BELOW 
DO  130  II-l.N 

130  P(II)*ALOG(X(II)-SIG) 

CALL  LIN(W,WM,P,EK,N,SIG,A,B,RES) 

DO  108  I-1,N 

PCAL-1 . 0-EXP ( -EC* ( X ( I ) -E  S IGO ) **EM ) 

108  WRITE(5 , 110)  PCAL,F( I) ,W(I) ,WM(I) 

110  F0RHAT(1X,°PCAL,F( I) ,W( I) ,WM(I)*° ,E12.6,3(1X,E12.6)) 

STOP 

END 

SUBROUTINE  LIN(W,WN,P,EK,N,SIG,A,B,RES) 

DIMENSION  W(20),P(20),EK(20),WN(20) 

C-0. 

D-0. 

E-0. 

G-0. 

H-0. 

DO  4 I-1,N 
C-C+(W(I)) 

D-D+(W(I)*P(I)) 

E-E+(W(I)*P(I)*P(I)) 

G«G+(W(I)*EK(I)) 

H*H+ ( W( I )*EK( I )*P ( I ) ) 

4 CONTINUE 
DEN-E*C-D*D 

A*( ( E*G-D*H) /( DEN ) ) 

B-((C*H-D*G)/(DEN)) 

C THE  CONSTANTS  A AND  B ARE  KNOWN 
C NOW  FIND  RESIDUE 
SUMM-0 . 

DO  5 1*1 ,N 

WN( I )■( A+B*P( I ) ) -EK( I ) 

5 SUMM*SUMM+W( I)*WN(I)*WN( I) 

RES-SUMM 

RETURN 

END 


C 

THIS 

IS  A BIAXIAL  PROGRAM 

00200 

IMPLICIT  DOUBLE  PRECIS  ION  (A-H-.O-Z ) 

00300 

DIMENSION  TH(20) , EM(20) ,W(20) ,A(20,20) ,B(20) ,XX(20) 

00400 

1 

,XK(20) 

00500 

c 

FOURIER  LEAST  SQUARE  APPROXIMATION  PROGRAM 

00600 

c 

ND-THE  FILE  NAME  AS  IN  ND-5(THIS  TERMINAL)  OR 

00700 

c 

FORND 

.DAT-THE  DISK  FILE 

00800 

WRITE(5 ,40) 

00900 

40 

FORMAT(2X, “ENTER  DATA  FILE  NO.0) 

01000 

RE AD (5 ,*)  ND 

01100 

WRITE(5 , 70) 

01200 

70 

FORMATC IX, “NOW  ENTER  NO.  OF  POINTS  °) 

01300 

READ(5 ,*)  N 

01400 

C 

N-NUMBER  OF  WEIBULL  POINT  PAIRS  OR  STRESS/FRACTURE/ 

01500 

C 

PERCENT  PAIRS. 

01600 

READ(ND,*)  ( (TH( I ) ,EM( I) ,W(I)) , 1-1 , N) 

01700 

DO  50  1*1  ,N 

01800 

50 

WRITE (5 ,60)  (I,TH(I),EM(I),W(I)) 

01900 

60 

FORMAT ( IX , “TH(  I)-,EM(  I)  , W( I)  FORI-“ , 12 ,3( 1X.E13 .6) ) 

02000 

C 

W-THE 

WEIGHT  AT  EACH  THETA 

02100 

C 

THETA 

IS  IN  RADIANS 

02200 

C 

EM-THE  WEIBUL  MODULUS 

02300 

C 

TO  MAKE  THE  MATRIX  SYMMETRIC  WE  MULTIPLY  TOP  ROW 

02400 

C 

BY  .5 

AFTER  FORMULATION  OF  MATRIX 

02500 

C 

DATAS 

ARE  UNEQUALLY  SPACED 

02600 

C 

FOURIER  APPROXIMATION  FITTING  THE  DATA  POINTS. 

02700 

C 

THEN  ' 

THE  NORMAL  EQUATIONS  ARE  SET  UP. 

02800 

C 

THEN  ' 

THEY  ARE  SOLVED  AND  TESTED  FOR  SIGNIFICANCE 

02900 

C 

AND  THE  FOURIER  APPROXIMATION  IS  INCREASED  BY  ONE  DEGREE 

03000 

e 

NEEDED  AND  THE  PROCESS  REPEATED  IF  NEEDED. 

03100 

c 

ITER-THE  NUMBER  OF  TIMES  THE  PROCESS  IS  EXECUTED. 

03200 

NP-5 

03300 

MP-N/2 

03400 

MP-2*MP 

03500 

ITER-0 

03600 

SIG0-30E+30 

03700 

M-l 

03800 

c 

SET  UP  THE  FIRST  MATRIX 

03900 

13 

1-1 

04000 

B( I)-PHO( I ,N,W, EM,TH) 

04100 

DO  51  I-2.M+1 

04200 

K-I-l 

04300 

51 

B(I)-RHO(K,N,W,EM,TH) 

04400 

DO  52  I-M+2 , 2*M+1 

04500 

K-I-M-l 

04600 

52 

B(I)-GHO(K,N,W,EM,TH) 

04700 

WRITE (5, 204) 

04800 

204 

FORMAT( / , 10X , “RIGHT  HAND  SIDE  MATRIX* ) 

04900 

WRITE(5 ,205)  (I,B(I) ,1-1 ,2*M+1) 

05000 

205 

FORMAT ( /,10X,I2I4X,E13.6) 

05100 

1-1 

05200 

J-l 

05300 

A(I,J)-RMC(I,J,N,W,EM,TH) 

05400 

1-1 

05500 

DO  53  J-1,M 

05600 

KO-I 

05700 

LO-J+1 

05800 

53 

A(K0,L0)-AMC(J,K0,L0,N,W,EM,TH) 

05900 

1-1 

06000 

DO  54  J-l.M 

06100 

K4-I 

06200 

L4-J+1+M 

06300 

54 

A(K4,L4)-BMC(J,K4,L4,N,W,EM,TH) 

06400 

DO  55  I-l.M 

06500 

DO  55  J-1,M 

06600 

Kl-I+1 

06700 

Ll-J+1 

06800 

55 

A(K1,L1)-DMC(I,J,K1,L1,N,W,EM,TH) 

06900 

DO  56  I-l.M 

07000 

DO  56  J-l.M 

07100 

K2-I+1 

07200 

L2-J+1+M 

07300 

56 

A( K2 , L2 ) -GMC ( I , J , K2 , L 2 , N , W , EM , TH ) 

07400 

DO  57  I-l.M 

07500 

DO  57  J-l.M 

07600 

K3-I+1+M 

07700 

L3-J+1+M 

07800 

57 

A(K3 ,L3)-HMC( I , J ,K3 ,L3 .N.W.EM.TH) 

07900 

DO  58  1-1 ,2*M 

08000 

DO  58  J-I+l , 2*M+1 

08100 

58 

A( J, I)-A( I , J) 

08200 

A(I,J)-A( J.I) 

08300 

WRITE(5 ,219) 

08400 

219 

FORMAT (/,6X, “MATRIX  OF  COEFICIENTS0/ ,7X,“I° ,5X, 

08500 

1 9X,“A(I, J)“ ,/) 

08600 

WRITE(5 , 220)  ( (I, J,  A(  I , J) , J-l , 2*M+1 ) , 1-1 , 2*M+1 ) 

08700 

220 

FORMAT(6X, I2.4X.I2, 5X.E13 .6) 

08800 

C NOW 

SOLVE  THE  NORMAL  MATRIX  EQUATIONS 

08900 

C NP 

IS  THE  NUMBER  OF  PRINTER 

09000 

M1«2*M+1 

09100 

CALL  S0LVE(A,XX,B,M1 .DET.NP) 

09200 

WRITE(5,225)  M,(XX(I) ,1-l.Ml) 

09300 

225 

FORMAT (2X,®M-° ,15,/ , lX,°XX(I)-° ,4(3E13.6)) 

09400 

C NOW 

TEST  THE  STATISTICAL  SIGNIFICANCE  OF  INCREASING  M 

09500 

DO  59  1-1 , 2*M+1 

09600 

59 

XK(I)-XX(I) 

09700 

IF(MP.EQ.N)  GO  TO  101 

09800 

IF(MP.LT.N)  GO  TO  102 

09900 

101 

IF(2*M+1 .EQ.(MP-l))  GO  TO  99 

10000 

102 

IF(2*M+1 . EQ.N)  GO  TO  99 

10100 

C COMPUTE  DELTA**2  THEN  SIGMA** 2 

10200 

DELN-0 . 

10300 

DO  61  1-1  ,N 

10400 

TEMP-XX(l)/2 

92 
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12700 
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DO  62  J-1,M 

Jl-J+1 

J2-M+1+J 

62  TEMP-TEMP+(XX(J1)*C0S(J*TH(I))+XX(J2)*SIN(J*TH(I))) 

61  DELN-DELN+W(I)*(EM(I)-TEMP)**2 

C NOW  COMPUTE  SIGMA** 2 
XT-N-2*M-1 
SIGN-DELN/XT 

WRITE (5, 200)  DELN.SIGN.SIGO 

200  FORMAT (2X , °DELN-DELTA**2-° , El 3 . 6/ , 

1 2X,oSIGN-SIGMA**2-°,E13.6,2X,#SIG0-°,E13.6) 

IF(SIGN-SIGO)  11,99,99 

11  SIGO-SIGN 

M-M+l 
GO  TO  13 

99  STOP 

END 

C FUNCTION  PHO  COMPUTES  THE  ELEMENTS  IN  R.H.S.  MATRIX 
FUNCTION  PHO(I,N,W,EM,TH) 

IMPLICIT  DOUBLE  PRECISION(A-H.O-Z) 

DIMENSION  W( 20 ) , EM( 20 ) , TH( 20 ) 

SUMM-0. 

DO  1 L-1,N 

1 SUMM-SUMM+(W(L)*EM(L))*.5 

PHO-SUMM 
RETURN 
END 

C FUNCTION  RHO  COMPUTES  THE  ELEMENTS  IN  R.H.S.  MATRIX 
FUNCTION  RHO ( J , N , W , EM , TH ) 

IMPLICIT  DOUBLE  PRECIS ION (A-H.O-Z) 

DIMENSION  W(20) ,EM(20) ,TH(20) 

SUMM-0. 

DO  1 L-l ,N 

1 SUMM-SUMM+ (W(L)*EM(L)*COS( J*TH ( L ) ) ) 

RHO-SUMM 

RETURN 

END 

C FUNCTION  GHO  COMPUTES  THE  ELEMENTS  IN  R.H.S.  MATRIX 
FUNCTION  GHO(J,N,W,EM,TH) 

IMPLICIT  DOUBLE  PRECISION(A-H ,0-Z) 

DIMENSION  W(20) ,EM(20) ,TH(20) 

SUMM-0. 


14700  DO  1 L-l ,N 

14800  1 SUMM-SUMM+(W(L)*EM(L)*SIN( J*TH(L) ) ) 

14900  GHO-SUMM 

15000  RETURN 

15100  END 

15200  C ALL  THE  FUNCTIONS  BELOW  FORMS  THE  L.H.S.  MATRIX 

15300  FUNCTION  RMC(I,J,N,W,EM,TH) 

15400  IMPLICIT  DOUBLE  PRECISION(A-H.O-Z) 

15500  DIMENSION  W(20) ,EM(20) ,TH(20) 

15600  SUMM-0. 


15700 

DO  1 L«1,N 

15800 

1 

SUMM-SUMM+( W(L) )* . 25 

15900 

RMC-SUMM 

16000 

RETURN 

16100 

END 

16200 

FUNCTION  AMC(J,KO,LO,N,W,EM,TH) 

16300 

IMPLICIT  DOUBLE  PRECIS ION (A-H.O-Z) 

16400 

DIMENSION  W(20) ,EM(20) ,TH(20) 

16500 

E-0. 

16600 

DO  1 L-1,N 

16700 

1 E-E+(W(L)*COS(J*TH(L)))*.5 

16800 

AMC-E 

16900 

RETURN 

17000 

END 

17100 

FUNCTION  BMC(J,K4,L4,N,W,EM,TH) 

IMPLICIT  DOUBLE  PRECISION (A-H.O-Z) 

17200 

17300 

DIMENSION  W(20),EM(20),TH(20) 

17400 

SUMM-0. 

17500 

DO  1 L-1,N 

17600 

1 

SUMM-SUMM+(W(L)*SIN(J*TH(L)))*.5 

17700 

BMC-SUMM 

17800 

RETURN 

17900 

END 

18000 

FUNCTION  DMC(I,J(K1,L1,N,W,EM,TH) 

IMPLICIT  DOUBLE  PRECISION(A-H.O-Z) 

18100 

18200 

DIMENSION  W(20) ,EM(20) ,TH(20) 

18300 

SUMM-0. 

18400 

DO  1 L-l ,N 

18500 

1 

SUMM-SUMM+ (W(L)*COS( J*TH(L) )*COS(I*TH(L) ) ) 

18600 

DMC-SUMM 

18700 

RETURN 

18800 

END 

18900 

FUNCTION  GMC(I, J,K2,L2,N,W,EM,TH) 

19000 

IMPLICIT  DOUBLE  PRECIS ION (A-H.O-Z) 

19100 

DIMENSION  W(20) , EM(20) ,TH(20) 

19200 

SUMM-0. 

19300 

19400 

DO  1 L-l ,N 

19500 

1 

SUMM-SUMM+ ( W(L )*S IN( J*TH( L ) )*COS ( I*TH( L ) ) ) 

19600 

GMC-SUMM 

19700 

RETURN 

19800 

END 

19900 

FUNCTION  HMC(I,J,K3,L3,N,W,EM,TH) 

20000 

IMPLICIT  DOUBLE  PRECISION(A-H.O-Z) 

20100 

DIMENSION  W(20),EM(20),TH(20) 

20200 

SUMM-0. 

20300 

DO  1 L-l , N 

20400 

1 

SUMM-SUMM+(W(L)*SIN( J*TH(L) )*SIN( I*TH(L) ) ) 

20500 

RMC-SUMM 

20600 

RETURN 

20700 

END 

94 


.00 

C 

PROGRAM  00046 ( INPUT , OUTPUT , TAPE5-INPUT , TAPE6-OUTPUT , TAPE 10 , TAPE 1 1 , 

00200 

C 

1TAPE 1 2 , TAPE 1 3 , TAPE 14 , TAPE 1 5 , TAPE 1 6 , PUNCH ) 

00300 

IMPLICIT  REAL*8 ( A-H , O-Z ) 

00400 

REAL*8  TITLE, WORD1 , WORD2 , WORD3 

00500 

REAL  NU21.NU31.NU32 

00600 

C**AMG054  FINITE  ELEMENT  STRESS  ANALYSIS  OF  ROCKET  NOZZLES 

00700 

INTEGER  CODE, ERR, PCODE , BW.TOPT 

00800 

DIMENSION  R(30 , 50) ,Z( 30 , 50) .CODE (30 , 50) ,TITLE( 13) ,UF(100),WF(100) 

00900 

1 , TANF( 100) , IMAX(50) ,IMIN(50) ,BC( 100 , 3) ,CONPR( 16 , 15) , IP(200) 

01000 

2 ,JP(200),P(200),TAU(200),PCODE(200),NEQ(50),  S(8,8) 

01100 

3 ,F(8) ,C(4 ,4) ,A( 120 ,80) ,B( 120) ,U(30 ,50) ,W(30,50) 

01200 

4 , RT(2000) ,ZT(2000) , TEMP (2000) ,PSI(4) ,SSMAX( 14) , SSMIN( 14) , IJSS( 14 

01300 

5 ,4) 

01400 

COMMON/TEM/A 

01500 

COMMON  UF,WF,TANF,JMIN,JMAX, ERR, MAX, R,Z, CODE, TITLE, IMAX.IMIN, 

01600 

1 CONPR, IP, JP,P,TAU, PCODE, BW,NEQ,JRAN,  S.F.NPCARD 

01700 

2 , TOPT,NP,MN 

01800 

COMMON/ VARPRO/ TABLE (12,10,5) 

01900 

EQUIVALENCE  (BC ,UF) , (A( 1 ) ,U) , (A( 1501 ) ,W) , (A( 3001 ), RT) 

02000 

1 , ( A(5001 ) ,ZT) , (A( 7001 ) ,TEMP) 

02100 

DATA  WORD1/6HH  ONLY/ .WORD2/6HEND  OF/ .WORD3/6HAMG054/ 

02200 

CALL  ERRSET( 209, 256, -1,1,0) 

02300 

KPLOT-O 

02400 

MAX  - 30 

02500 

JR1  - 120 

02600 

READ(40 , 101 ) N1.N2.N3 

02700 

C Nl- 

■INPUT  DATA  FILE  NO. -40 

02800 

C N2« 

■OUTPUT  DATA  FILE  NO. -41 

02900 

C N3- 

■MAT  NO  FOR  WHICH  WEIBULL  STS.  IS  TO  BE  DONE! FOR  WEIBUL 

03000 

4 

■ READ(40 , 100 )(TITLE( I) ,1-1,13) 

03100 

IF( TITLE ( 1 3 ) . NE . WORD1 ) WRITE (46 ) TITLE 

03200 

100C 

1 READ(40 , 101 ) JMIN.JMAX,  NCONT , TOPT , NTABLE 

03300 

REWIND  42 

03400 

REWIND  44 

03500 

REWIND  45 

03600 

CALL  ERRSET(l) 

03700 

WRITE(41 ,205)  (TITLE(I) ,1-1 , 13) 

03800 

WRITE (4 1,206)  JMIN.JMAX,  NCONT, NTABLE 

03900 

ITOPT-TOPT+1 

04000 

GO  TO( 3000 , 3001 , 3002 , 3003 , 3004 , 3000 ) , ITOPT 

04100 

3001 

WRITE(41 , 218) 

04200 

GO  TO  3000 

04300 

3002 

WRITE(41 , 21 9) 

04400 

GO  TO  3000 

04500 

3003 

WRITE(41 , 220) 

04600 

GO  TO  3000 

04700 

3004 

WRITE(41 , 221 ) 

04800 

3000 

CONTINUE 

04900 

JRAN-JMAX-JMIN+1 

05000 

c*****  initialize 

05100 

DO  2 J-l.JRAN 

05200 

DO  1 1-1 , MAX 
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05300 

R(I,J)«0. 

05400 

Z(I,J)-0. 

05500 

1 

CODE(I,J)«0 

05600 

IMIN(J)-  1000 

05700 

2 

IMAX(J)-  0 

05800 

DO  3 1-1,100 

05900 

DO  3 J-1,3 

06000 

3 

BC(I, J)-0. 

06100 

ERR-0 

06200 

c*** 

READ  NODAL  POINT  DATA 

06300 

CALL  MESH 

06400 

c 

CALL  GRIDSC( R , Z , CODE , IMIN , IMAX , JRAN , KPLOT , TITLE ) 

06500 

IF(TOPT.EQ.5)  TOPT-O 

06600 

IF(ERR.EQ.O)  GO  TO  6 

06700 

WRITE(41 ,212)  ERR 

06800 

5 

READ(40 , 100)  (TITLE(I), 1-1,13) 

06900 

IF(WORD2.NE.TITLE(l).AND.WORD3.NE.TITLE(l))  GO  TO 

07000 

GO  TO  2001 

07100 

6 

IF  (TITLE (13)  .EQ.  WORDl)  GO  TO  5 

07200 

c*** 

READ  PRESSURE  CARDS 

07300 

CALL  PRESBC( IP , JP , P , TAU , PCODE , NPCARD , ERR) 

07400 

IF(ERR.NE.O)GO  TO  5 

07500 

IF ( NPCARD. EQ.O)  GO  TO  8 

07600 

WRITE(41 , 207 ) 

07700 

c*** 

ORDER  THE  PRESSURE  CARDS  BY  INCREASING  I AND  J 

07800 

IF( NPCARD. EQ.l)  GO  TO  706 

07900 

N1 -NPCARD- 1 

08000 

DO  704  N-l.Nl 

08100 

IS-  IP(N) 

08200 

JS-JP(N) 

08300 

NS-N 

08400 

Ml-  N+l 

08500 

DO  703  M-Ml, NPCARD 

08600 

IF(JP(M)-JS)702,701 ,703 

08700 

701 

IF(IP(M).GT.IS)GO  TO  703 

08800 

702 

IS-IP(M) 

08900 

JS-JP(M) 

09000 

NS-M 

09100 

703 

CONTINUE 

09200 

TS-TAU(NS) 

09300 

PS-P(NS) 

09400 

IPCS-PCODE(NS) 

09500 

IP(NS)-IP(N) 

09600 

JP(NS)-JP(N) 

09700 

P(NS)-P(N) 

09800 

TAU(NS)-TAU(N) 

09900 

PCODE (NS )-PCODE(N) 

10000 

IP(N)-IS 

10100 

JP(N)-JS 

10200 

P(N)-PS 

10300 

TAU(N)-TS 

10400 

PCODE (N)-IPCS 
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10500 

704  WRITE (41 , 208)IS, JS,PS,TS, IPCS 

10600 

706 

WRITE(41 ,208)IP(NPCARD) , JP(NPCARD) ,P(NPCARD) .TAU(NPCARD) ,PCODE(NPCA 

10700 

1 RD) 

10800 

c*** 

READ  CONTINUUM  MATERIAL  PROPERTIES 

10900 

8 

IF(NCONT.EQ.O)  GO  TO  13 

11000 

WRITE(41 ,203) 

11100 

DO  9 N-l.NCONT 

11200 

READ(40,103)  MN.El ,E2,E3,NU21 ,NU31 ,NU32,PHI ,G13,ALPHA1 , 

11300 

1 ALPHA2 , ALPHA3 , CONPR ( 1 5 , MN ) , CONPR ( 1 6 , MN ) 

11400 

WRITE (41 ,204)  MN,E1 ,NU21 , ALPHA1 ,G13 ,PHI ,E2,NU31 , ALPHA2 , 

11500 

i 

1 CONPR(15,MN),E3,NU32,ALPHA3,CONPR(16,MN) 

11600 

CALL  PR0P(E1 ,E2 ,E3,NU21 ,NU31 ,NU32,G13 ,PHI , ALPHA1 , 

11700 

1 ALPHA2,ALPHA3,C,PSI) 

11800 

ICOUNT-O 

11900 

DO  801  11-1,4 

12000 

C0NPR(II+10,MN)-PSI(II) 

12100 

DO  801  JJ-11,4 

12200 

ICOUNT-ICOUNT+1 

12300 

801 

CONPR( ICOUNT ,MN)-C( II , JJ ) 

12400 

9 

CONTINUE 

12500 

13 

IF(NTABLE.EQ.O)GO  TO  11 

12600 

c*** 

READ  MATERIAL  PROPERTY  VS.  TEMPERATURE  TABLES 

12700 

WRITE(41 ,215) 

12800 

DO  10  N-l , NT ABLE 

12900 

READ(40, 104)  MTN , NTEMP , BFR , BFZ , PHI , TO 

13000 

WRITE (41 ,216)  MTN, BFR, BFZ, PHI 

13100 

MTN  - MTN  - 15 

13200 

DO  1001  Nl-1, NTEMP 

13300 

RE AD(40,105)(TABLE(II,N1, MTN), 11-1,11) 

13400 

1001 

WRITE(41 ,217)  (TABLE(II,N1, MTN), 11-1,11) 

13500 

TABLE ( 1 2 , 1 , MTN ) -NTEMP 

13600 

TABLE (12,2, MTN ) -BFR 

13700 

TABLE (12, 3, MTN) -BFZ 

13800 

TABLE(12 ,4,MTN)»MTN+  15 

13900 

TABLE( 12 , 5 ,MTN)»T0 

14000 

TABLE(12,6,MTN)  - PHI 

14100 

10 

CONTINUE 

14200 

c*** 

DETERMINE  TEMPERATURES  FOR  CONTINUUM  ELEMENTS 

14300 

11 

REWIND  42 

14400 

IF(TOPT.NE.O)CALL  TEMPT ( TO PT) 

14500 

REWIND  42 

14600 

CALL  SETUP(A,B) 

14700 

c*** 

SOLVE  FOR  DISPLACEMENTS 

14800 

CALL  BACSUB(A,B,NEQ,BW, JRAN) 

14900 

REWIND  45 

15000 

DO  14  J-l.JRAN 

15100 

J 1-JR1  -JRAN+J 

15200 

NEQJ  - NEQ(J) 

15300 

14 

WRITE (45)  (A(J1,N1), Nl-1, NEQJ) 

15400 

REWIND  45 

15500 

WRITE(41 ,213) 

15600 

DO  15  J-l.JRAN 

97 


15700 

N2-  NEQ( J)/2 

15800 

READ(45)  (U(N1,J),W(N1 

,J)  .N1-1.N2) 

15900 

11-  IMIN(J) 

16000 

12-  IMAX(J) 

16100 

Jl-  JMIN-l+J 

16200 

DO  15  Z-11,12 

16300 

13-  I+1-IMIN(J) 

16400 

NNN-0 

16500 

IF(  J . EQ . JRAN . AND . I . EQ . 1 2 )NNN—  1 

16600 

WRITE(46)  R(I3,J),Z(I3,J),U(I3,J),W(I3,J),NNN 

16700 

RPDR-R(I3,J)+U(I3,J) 

16800 

ZPDZ-Z(I3,J)+W(I3,J) 

16900 

15  WRITE(41 ,214)  I,J1,  R(I3,J),  Z(I3,J),  U(I3,J),  W(I3,J),  RPDR,  ZPDZ 

17000 

c****  calculation  OF  ELEMENT 

STRESSES 

17100 

REWIND  44 

17200 

JRAN1-  JMAX-JMIN 

17300 

IC-0 

17400 

DO  1501  1-1,14 

17500 

SSMAX(I)— 1.E20 

17600 

1501  SSMIN(I)-1.E20 

17700 

DO  16  J-1.JRAN1 

17800 

Il-IMIN(J) 

17900 

I2-IMAX(J)-1 

18000 

• DO  16  1-11,12 

18100 

I3-I+1-I1 

18200 

J3-JMIN-1+J 

18300 

IN-I-IMIN(J+1)+1 

18400 

16  IF( CODE ( 13 , J ) / 1000000 . LE . 25 ) CALL  STRESS(U,W,R,Z,I3,J,I,J3,CODE,IC, 

18500 

1 IN , SSMAX , SSMIN , I JSS ) 

18600 

XX-1.E20 

18700 

NNN— 1 

18800 

WRITE (46 ) 13 , JRAN ,( XX , I- 

1,16), NNN 

18900 

WRITE (41 , 222 ) ( IJSS( I , 1 ) , IJSS( I , 2) , SSMIN( I ) , I JSS( I , 3) , 

19000 

1 I JSS (I, 4) , SSMAX(I) ,1 

-1,14) 

19100 

222  FORMAT ( 1 H 1 1 2X5 9HMINIMUM  AND  MAXIMUM  VALUES  OF  STRESS  AND  STRAIN  IN 

19200 

1 THE  BODY/1H041X7HMINIMUM23X , 7HMAXIMUM/ 1H0 , 8HQUANTITY25X1HI , 4X1HJ 

19300 

21 0X5HVALUE 14X1HI , 4X1HJ 10X5HVALUE/ 

19400 

317H0RADIAL  STRESS 

13X215 , 1PE15 .4, 10X215 ,E15 .4/ 

19500 

417HOHOOP  STRESS 

13X215, E15. 4, 10X215, E15. 4/ 

19600 

517HOAXIAL  STRESS 

13X215, E15. 4, 10X215, E15. 4/ 

19700 

617H0R-Z  SHEAR  STRESS 

13X215, E15. 4, 10X215, E15. 4/ 

19800 

717H0MAX  STRESS 

13X215, E15. 4, 10X215, E15. 4/ 

19900 

817H0MIN  STRESS 

13X215, E15. 4, 10X215, E15. 4/ 

20000 

917H0MAX  SHEAR  STRESS 

13X215 ,E15 .4,10X215 ,E15 .4/ 

20100 

117H0RADIAL  STRAIN 

13X215, E15. 4, 10X215, E15. 4/ 

20200 

1 1 7HOHOOP  STRAIN 

13X215, E15. 4, 10X215, E15. 4/ 

20300 

217HOAXIAL  STRAIN 

13X215, E15. 4, 10X215, E15. 4/ 

20400 

317HOR-Z  SHEAR  STRAIN 

13X215, E15. 4, 10X215, E15. 4/ 

20500 

417H0MAX  STRAIN 

13X215, E15. 4, 10X215, E15. 4/ 

20600 

517H0MIN  STRAIN 

13X215, E15. 4, 10X215, E15. 4/ 

20700 

617H0MAX  SHEAR  STRAIN 

13X215, E15. 4, 10X215, E15. 4) 

20800 

READ(40, 100)  (TITLE(I) , 

1-1,13) 
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20900 

21000 

21100 

21200 

21300 

21400 

21500 

21600 

21700 

21800 

21900 

22000 

22100 

22200 

22300 

22400 

22500 

22600 

22700 

22800 

22900 

23000 

23100 

23200 

23300 

23400 

23500 

23600 

23700 

23800 

23900 

24000 

24100 

24200 

24300 

24400 

24500 

24600 

24700 

24800 

24900 

25000 

25100 

25200 

25300 

25400 

25500 

25600 

25700 

25800 

25900 

26000 


2001  CONTINUE 

IF(W0RD1 . NE . TITLE (13)) WRITE ( 46 )TITLE 
IF(WORD2.NE.TITLE(l))  GO  TO  1000 
CALL  CLEAN 
CALL  EXIT 

100  FORMATU3A6) 

101  FORMAT(8I5) 

102  FORMAT( 215 , 2F10 . 5 , 15 ) 

103  FORMAT(I10,6F10.5/7F10.5) 

104  FORMAT(2I10,4F10.5) 

105  F0RMAT(7F10.5/10X4F10.5) 

201  FORMAT(2I10 , 2F20.4, 110) 

203  F0RMATQH125X38HM  ATERIAL  PROPERTIES) 

204  FORMAT (10H0MAT.  NO.-I2.4X3HE1-1PE11 .4.4X5HNU21-E11 .4, 2X7HALPHA1- 

1 Ell .4.5X4HG13-E11 .4.5X4HPHI-E11 .4/ 

2 16X3HE2-E11 .4.4X5HNU31-E11 .4.2X7HALPHA2" 

2 E 1 1 . 4 , 5X4HBFR-E 1 1 . 4/ 1 6X3HE  3-E 1 1 . 4 , 4X5HNU3  2-E 1 1 . 4 , 2X7HALPHA3- 

3 E11.4.5X4HBFZ-  Ell. 4) 

205  FORMAT(49H1AMG054  FINITE  ELEMENT  STRESS  ANALYSIS  OF  NOZZLES/ 

184H0WRITTEN  BY  JOHN  BRISBANE,  ROHM  AND  HAAS  CO.,  REDSTONE  ARSENAL 
2 RESEARCH  LABORATORIES/ 1H013A6) 

206  FORMAT ( 6H0JMIN-I5/6H  JMAX-I5/ 

1 18H  NO.  OF  MATERIALS-I5/ 

225H  NO.  OF  MATERIAL  TABLES  - 13) 

207  F0RMAT(1H1 16X34 HA  PPLIED  PRESSURES/ 
19X1HI9X1HJ13X1HP17X3HTAU11X5HPCODE  ) 

208  F0RMAT(2I10, 2F20 .5 , 110) 

212  FORMAT(I5,52H  DATA  ERRORS  NOTED,  PROGRAM  PROCEEDS  TO  NEXT  PROBLEM) 

213  FORMAT ( 1H13X1HI4X1HJ6X12HR-COORDINATE6X12HZ-COORDINATE6X14HR-DISPL 
1ACEMENT6X14HZ-DISPLACEMENT6X1 1HR  ♦ DELTA-R6X1 1HZ  + DELTA-Z  /) 

214  FORMAT ( 215,  2F15.4,  2(  10X,  1PE15.4),  0P2F15.4) 

215  FORMAT ( 1H1,35X,49HM  ATERIAL  PROPERTY  TABLE 
l S) 

216  F0RMAT(1H0/13H  MATERIAL  NO. 13, 10X4HBFR-1PE13.6, 10X4HBFZ-E13 .6, 

1 10X4HPHI«E13.6/1H0,7X,1HT,7X,2HE1,10X,2HE2,10X,2HE3,8X,4RNU21, 

2 4X , 4HNU3 1 , 4X , 4HNU32 , 6X , 3HG1 3 , 8X , 6HALPHA1 , 6X , 6HALPHA2 , 6X , 6HALPHA3 ) 

217  F0RMAT(1H  F9.0,3E12.4,3F8.5,4E12.4) 

218  FORMAT (35H0THE  BODY  HAS  A UNIFORM  TEMPERATURE) 

219  FORMATC40H0THE  TEMPERATURE  TABLE  WAS  INPUT  ON  TAPE) 

220  FORMAT (4 1H0THE  TEMPERATURE  TABLE  WAS  INPUT  ON  CARDS) 

221  FORMAT ( 6 8HONODAL  POINT  TEMPERATURES  WERE  INPUT  ON  TAPE  AS  DETERMIN 
1ED  BY  AMG065) 

STOP 

END 

SUBROUTINE  PR0P(E1 , E2 , E3 , NU21 , NU31 , NU32 , G13 , PHI , ALPHAl , 

1 ALPHA2,ALPHA3,C,PSI) 

IMPLICIT  REAL*8 ( A-H , O-Z ) 

REAL*8  TITLE, WORD1 .WORD2.WORD3 
REAL  NU21.NU31.NU32 

DIMENSION  C(4,4) ,PSI(4) ,CD(4,4) ,D(4,4) 

DATA  D(2,2) ,D(1 ,2) ,D(3,2) ,D(4, 2) ,D(2,1),D(2,3),D(2,4)/1. ,6*0./ 

C***  FORM  THE  C MATRIX 


sin(x)-dsin(x) 

COS(X)-DCOS(X) 

CFU  - l./E2/E3-NU32*NU32/E2/E2 

CF12  - -NU21/E1/E3-NU31*NU32/EI/E2 

CF13  - NU21*NU32/E1/E2+NU31/E1/E2 

CF22  - 1./E1/E3-NU31*NU31/E1/E1 

CF23  - -NU32/E1/E2-NU21*NU31/E1/E1 

CF33  - 1./E1/E2-NU21*NU21/E1/E1 

DET  - 1./E1*CF1H-NU21/E1*CF12-NU31/E1*CF13 

C(1,1)-CF11/DET 

C(l,2)— CF12/DET 

C(1,3)-CF13/DET 

C(2,2)-CF22/DET 

C(2,3)— CF23/DET 

C(3,3)-CF33/DET 

C(4,4)-C13 

C(l,4)-0. 

C(2,4)-0. 

C(3,4)-0. 

C(2,l)-C(l,2) 

C(3,l)-C(l,3) 

C(4,l)»0. 

C(3,2)-C(2,3) 

C(4,2)»0. 

C(4,3)-0. 

C***  FORM  THE  VECTOR  PS I 
DO  1 N-1,3 

1 PSI(N)-C(N,1)*ALPHA1+C(N,2)*ALPHA2+C(N,3)*ALPHA3 
PSI(4)-0. 

IF(PHI .EQ.O. ) GO  TO  8 
C***  ROTATE  C AND  PS I BY  THE  ANGLE  PHI 
PHI-PHI/57.29578 
CP-COS ( PHI ) 

SP-SIN(PHI) 

D(1,1)-CP*CP 

D(1,3)-SP*SP 

D(1,4)-SP*CP 

D(3,l)-D(l,3) 

D(3,3)-D(l ,1) 

D(3,4)—  D(l,4) 

D(4,1)-2.*D(3,4) 

D(4,3)— D(4,l) 

D(4,4)-D(l ,1)-D(1 ,3) 

DO  2 1-1,4 
CD(I,l)-0. 

DO  2 K-1,4 

2 CD(I,1)-CD(I,1)+D(K,I)*PSI(K) 

DO  3 1-1,4 

3 PSI(I)-CD(I,1) 

DO  5 1-1,4 

DO  5 J-1,4 
CD(I, J)-0. 


100 


31300 

31400 

31500 

31600 

31700 

31800 

31900 

32000 

32100 

32200 

32300 

32400 

32500 

32600 

32700 

32800 

32900 

33000 

33100 

33200 

33300 

33400 

33500 

33600 

33700 

33800 

33900 

34000 

34100 

34200 

34300 

34400 

34500 

34600 

34700 

34800 

34900 

35000 

35100 

35200 

35300 

35400 

35500 

35600 

35700 

35800 

35900 

36000 

36100 

36200 

36300 

36400 


DO  5 K-1,4 

5 CD(I,J)-CD(I,J)+C(I,K)*D(K,J) 

DO  7 >1,4 

DO  7 J-I.4 
C(I,J)-0. 

DO  6 K-1,4 

6 C(I, J)-C(l, J)+D(K, I)*CD(K, J) 

7 C(J,I)«C(I,J) 

8 RETURN 
END 

SUBROUTINE  STIFFQ(I,J,PI,PT,IPT) 

C*****  CALCULATION  OF  STIFFNESS  MATRIX  OF  A QUADRILATERAL  COMPOSED 
C OF  FOUR  TRIANGLES. 

C CENTER  NODAL  DISPLACEMENTS  AND  ELEMENT  PRESSURE  ARE  ELIMINATED 

C AND  DATA  FOR  THEIR  CALCULATION  IS  WRITTEN  ON  TAPE  12 

IMPLICIT  REAL*8 ( A-H , O-Z ) 

REAL*8  TITLE , W0RD1 ,WORD2 , WORD3 
INTEGER  ERR, CODE,  TOPT 

COMMON  UF,WF,TANF,JMIN,JMAX, ERR, MAX, R,Z, CODE, TITLE, IMAX.IMIN, 

1 CONPR(16, 15) ,DUM1 ,S,F,NPCARD,TOPT,NP,MN 

DIMENSION  UF( 100) ,WF( 100) ,TANF( 100) ,R(30,50) ,Z( 30,50) .CODE (30, 50) , 

1 TITLE(13) ,IMAX(50) ,IMIN(50) ,DUM1(  726),  S(8,8),F(8) 

2 , IT(6 ,4) ,AK(10, 10) ,F2(2) ,RR(5 ) ,C(4,4) ,ET(4) ,ZZ(5),F1(10) 

DOUBLE  PRECISION  AK,F1,F2,DET 

DATA  IT(1 , 1 ) , IT(2 , 1 ) ,IT(3 , 1 ) ,IT(4, l) ,IT(5 , 1) ,IT(6 , 1) , IT( 1 , 2) , 

1 IT (2,2) , IT(3 ,2) , IT(4,2) , IT(5 , 2) , IT(6 ,2) ,IT(1 ,3) ,IT(2,3) ,IT(3 ,3) , 
2IT(4 , 3) , IT(5 ,3) , IT(6, 3) , IT( 1,4), IT(2 ,4) , IT(3 ,4) ,IT(4,4) ,IT(5 ,4) , 
3IT(6,4)/1,2,3,4,9,10,3,4,7,8,9,10,7,8,5,6,9, 10,5,6,1,2,9,10/ 

1000  FORMAT( 7H0STIFFQ) 

II-  I+IMIN(J)-1 

IN-I1-IMIN(J+1)+1 

RR(1)-R(I,J) 

RR(2)-R(I+1 , J) 

RR(3)-R(IN+1 , J+l) 

RR(4)-R( IN , J+l ) 

ZZ(1)-Z(I,J) 

ZZ(2)-Z(>1,J) 

ZZ(3)-Z(IN-*-l,J+l) 

ZZ(4)»Z(IN,J+1) 

RK-  (RR(l)+RR(2)+RR(3)+RR(4))/4.0 
ZK-  (ZZ(l)+ZZ(2)+ZZ(3)+ZZ(4))/4.0 
RR(5)-RR(1) 

ZZ(5)-ZZ(1) 

DT-0. 

IF( TOPT. NE.O) READ ( 42 )DT 

1001  PORMAT(8HOSTIFFQ1 ,110) 

IF(MN.LT. 16)  GO  TO  9 

CALL  INTERP(C,ET,BFR,BFZ,DT,MN) 

1002  FORMAK7HO INTER?) 

ET(1)  - ET(1)*DT 
ET(2)  - ET(2)*DT 
ET(3)  - ET(3)*DT 


101 


36500 

ET(4)  - ET(4)*DT 

36600 

GO  TO  10 

36700 

9 

ICOUNT-O 

36800 

DO  8 11-1,4 

36900 

ET(II)-CONPR(II+10,MN)*DT 

37000 

DO  8 JJ-11,4 

37100 

ICOUNT-ICOUNT+1 

37200 

C( II , JJ )-CONPR( 1COUNT ,MN) 

37300 

8 

C(JJ,II)-C(II,JJ) 

37400 

BFR-CONPR(15,MN) 

37500 

BFZ-CONPR(16,MN) 

37600 

10 

IPT1-IPT/10 

37700 

IPT2-IPT-10*IPT1 

37800 

DO  6 N-1,10 

37900 

DO  7 M-1,10 

38000 

7 

AK(N,M)-0. 

38100 

6 

Fl(H)-0. 

38200 

DO  1 N-1,4 

38300 

P-0. 

38400 

TAU-O. 

38500 

RI-RR(N) 

38600 

ZI-ZZ(H) 

38700 

RJ-RR(N+1) 

38800 

ZJ-ZZ(N+1) 

38900 

IF( IPT1 .EQ.N.OR. IPT2.EQ.N)  P-PI 

39000 

IF(IPT1 .EQ.N.0R.IPT2.EQ.N)  TAU-PT 

39100 

1003 

FORMAT*  9HOGOSTIFF3 ) 

39200 

13 

CALL  STIFF3(RI,RJ,RK,ZI,ZJ,ZK,P,TAU,IPT,C,ET,BFR,BFZ) 

39300 

1004 

FORMAT( 10H0REST1FF3  ) 

39400 

DO  1 Ml-1,6 

39500 

H1-IT(M1,N) 

39600 

F1(M1)-F1(N1)+F  (Ml) 

39700 

DO  1 M2-1.6 

39800 

N2-IT(M2,N) 

39900 

1 

AK(N1 ,N2)-AK(N1 , N2)+S  (Ml, M2) 

40000 

DET-AK( 9 , 9 )*AK( 10 , 10)-AK(9 , 10)**2 

40100 

1005 

FORMAT ( 5H0DET- , E20 . 8 ) 

40200 

AJC(9, 10)-AK(9,9)/DET 

40300 

AK(9,9)-AK(10,10)/DET 

40400 

AK(10,10)-AK(9,10) 

40500 

AK(10,9)— AK(10,9)/DET 

40600 

AK(9,10)-AK(10,9) 

40700 

DO  3 R-9,10 

40800 

F2(M-8)-0. 

40900 

DO  2 M-1,8 

41000 

AK(M,M)-0. 

41100 

DO  2 Hl-9,10 

41200 

2 

AK(H,M)-AK(N,M)+AK(N,N1)*AK(M,N1) 

41300 

DO  3 M-1,2 

41400 

3 

F2(H-8)-F2(H-8)+AR(N,M+8)*Fl(Mi-8) 

41500 

DO  5 N-1,8 

41600 

DO  401  M-1,8 
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41700 

DET-AK(N.M) 

41800 

DO  4 Nl-9,10 

41900 

* 

4 DET-DET-AK(N,N1)*AK(N1,M) 

42000 

401  S(N,M)-DET 

42100 

DO  402  M-1,2 

42200 

402  F1(N)-F1(N)-AK(N,M+8)*F2(M) 

42300 

5 F(N)-F1(N) 

42400 

WRITE (44)  C,ET,RK,ZK,F2,((AK(I1,I2),I1-9,10),I2-1,8) 

42500 

PI-0. 

42600 

PT-0. 

42700 

IPT-0 

42800 

1006  FORMAT(9HORESTIFFQ) 

42900 

RETURN 

43000 

END 

43100 

SUBROUTINE  INTERP( C , ET , BFR , BFZ , TEH  ,MN) 

43200 

c*** 

INTERPERTING  IN  MATERIAL  PROPERTY  TABLE 

43300 

IMPLICIT  REAL*8(A-H,0-Z) 

43400 

REAL*8  TITLE , WORD1 , WORD2 , WORD3 

43500 

DIMENSION  C(4,4),ET(4) , TABLE ( 12, 10) ,X(10) 

43600 

COMMON/VARPRO/XTABLE( 1 2 , 10,5) 

43700 

EQUIVALENCE  (X,E1 ) , (X(2) ,E2) , (X(3) ,E3) , (X(4) ,XNU21 ) ,(X(5 ) ,XNU31 ) , ( 

43800 

1 X(6) ,XNU32) ,(X(7),G13),(X(8) ,ALPHA1 ) ,(X(9) , ALPHA2) ,(X( 10) , ALPHA3 

43900 

2 ) 

44000 

DO  8 11-1,12 

44100 

^ DO  8 JJ-1,10 

44200 

8 

TABLE ( I I , JJ ) -XTABLE ( I I , J J , MN- 1 5 ) 

44300 

1 

NTEMP-TABLE (12 , 1 )+0 . 5 

44400 

TEMP-TEM  +TABLE (12,5) 

44500 

IF ( TEMP. LE .TABLE (1 , l))GO  TO  5 

44600 

IF( TEMP. GE. TABLE (l,NTEMP))GO  TO  6 

44700 

DO  2 N-2.NTEMP 

44800 

Nl-N 

44900 

IF ( TEMP. LE. TABLE (1 ,N) . AND. TEMP. GE. TABLE (1 ,N-l))GO  TO  3 

45000 

2 

CONTINUE 

45100 

3 

RATIO- (TEMP-TABLE ( 1 ,N1-1 ) ) / (TABLE( 1 ,N1 )-TABLE( 1 ,N1-1 ) ) 

45200 

7 

CONTINUE 

45300 

DO  4 11-2,11 

45400 

4 

X( II-l )-TABLE ( I I , N1 -1 ) +RATIO*( TABLE ( I I , N1 J-TABLE ( I I , N1 -1 ) ) 

45500 

PHI-TABLE (12, 6) 

45600 

BFR-TABLE(12,2) 

45700 

BFZ-TABLE(12,3) 

45800 

CALL  PR0P(E1 ,E2,E3,XNU21 ,XNU31 ,XNU32,G13,PHI,ALPHA1 , ALPHA2 , ALPHA3 , 

45900 

1 C,ET) 

46000 

RETURN 

46100 

5 

Nl-2 

46200 

RATIO-O . 

46300 

GO  TO  7 

46400 

6 

Nl-NTEMP 

46500 

RATIO-1 .0 

46600 

GO  TO  7 

46700 

END 

46800 

SUBROUTINE  MESH 
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D 


I 


46900 

47000 

47100 

47200 

47300 

47400 

47500 

47600 

47700 

47800 

47900 

48000 

48100 

48200 

48300 

48400 

48500 

48600 

48700 

48800 

48900 

49000 

49100 

49200 

49300 

49400 

49500 

49600 

49700 

49800 

49900 

50000 

50100 

50200 

50300 

50400 

50500 

50600 

50700 

50800 

50900 

51000 

51100 

51200 

51300 

51400 

51500 

51600 

51700 

51800 

51900 

52000 


IMPLICIT  REAL*8(A-H,0-Z) 

REAL*4  R1,Z1,W1,W2,W3 
REAL*4  POLAR,  SAME 
REAL*8  TITLE , WORDl , W0RD2 , W0RD3 
REAL  LINE 

INTEGER  CT, CODE, TYPE, ERR, PCODE.BW 
DIMENSION  Rl(30) ,Z1(30) 

DIMENSION  R(30,50) ,Z(30 ,50) , CODE (30, 50) ,TITLE( 13) ,UF(100),WF(100) 

1 , TANF( 100) , IMAX(50) ,IMIN(50) ,BC( 100 ,3) ,CONPR(16, 15) , IP (200) 

2 ,JP(200),P(200),TAU(200),PCODE(200),NEQ(50),  S(8,8) 

3 ,F(8),RR(5),ZZ(5) 

COMMON  UP, WF,TANF,JMIN,JMAX, ERR, MAX, R,Z, CODE, TITLE, IMAX.IMIN, 

1 CONPR , IP , JP , P , TAU , PCODE , BW , NEQ , JRAN , S.F.NPCARD 

2 , ITOPT,NP,MN 
EQUIVALENCE  (BC(1),UF) 

DATA  POLAR, LINE, SAME/1HP, 1HL, 1HS/ 

DATA  WORD1/6HH  ONLY/ 

C****  READ  NODAL  POINT  DATA 
ABS(X)-DABS(X) 

SIN(X)-DSIN(X) 

COS(X)-DCOS(X) 

WRITE (4 1,204) 

ICOUNT-O 

IC-1 

IF(IT0PT.NE.4.AND.  ITOPT  .NE.  5)  GO  TO  1 
JRAN-JMAX+l-JMIN 

READ(42)  (XMAX(J) ,IMIN( J) ,NEQ( J) , J-l , JRAN) ,BW 

BW-2*BW 

DO  1000  J-l, JRAN 
IRAN-NEQ(J) 

NEQ(J)  - 2*NEQ( J ) 

READ(42 )( R1 ( I ) ,Z1 (I) , C0DE(I , J) , 1-1 , IRAN) 

DO  1100  1-1 , IRAN 
R(I,J)-R1(I) 

1100  Z(I,J)-Z1(I) 

1000  CONTINUE 

1001  READ(40, 104)  W2,I, J.TYPE, II , 12, 13, I4,BC1 ,BC2,BC3,BC4 
IF(I.EQ.O)  GO  TO  9 

IF(J.LE.JMAX  ) GO  TO  1002 

ERR-ERR+1 

WRITE (41 , 200)  I,J 

1002  CT-I4+10*I3+100*I2+1000*I1 

IF(I1+I2+I3+I4.NE .0)  WRITE(41 ,203)  I,J,I1,BC1 ,I2,BC2,I4,BC4 

IF(ABS(BC1)+ABS(BC2)+ABS(BC3)+ABS(BC4) . EQ.O. ) GO  TO  1003 

UF(IC)-BC1 

WF( IC)-BC2 

TANF(IC)-BC4 

CT-CT+IC*10000 

IC-IC+1 

IF(IC.LE.IOO)  GO  TO  1003 

ERR-ERR+1 

WRITE(41 ,201 ) 


j 


j 


r 
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52100 

52200 

52300 

52400 

52500 

52600 

52700 

52800 

52900 

53000 

53100 

53200 

53300 

53400 

53500 

53600 

53700 

53800 

53900 

54000 

54100 

54200 

54300 

54400 

54500 

54600 

54700 

54800 

54900 

55000 

55100 

55200 

55300 

55400 

55500 

55600 

55700 

55800 

55900 

56000 

56100 

56200 

56300 

56400 

56500 

56600 

56700 

56800 

56900 

57000 

57100 

57200 


GO  TO  19 

1003  CT-CT+1000000*TYPE 
Jl-J-JMIN+1 
I8-I-IMIN(J1)+1 

CODE ( 18 , J1 ) - ( CODE ( 18 , J 1 ) / 1 00000000 )* 1 00000000+CT 
IF(W2.NE.LINE)  GO  TO  1001 
READ(40 , 104)  W2,II,JJ 
NSTEPS«MAXO(IABS( II-I ) , IABS( JJ-J ) ) 

ISTEP-( II-I ) /NSTEPS 
JSTEP-(JJ-J)/NSTEPS 
DO  1004  N-l, NSTEPS 
I-I+ISTEP 
J-J+JSTEP 

IF(I1+I2+I3+I4.NE.0)  WRITE (41 , 203 ) I , J, II ,BC1 , I2.BC2, I4,BC4 

Jl-J+l-JMIN 

I8-I-IMIN(J1 )+l 

1004  CODE ( 18 , J 1 ) - ( CODE ( 18 , J 1 ) / 1 00000000 )* 100000000+CT 
GO  TO  1001 

104  FORMAT ( IX, A1 ,13,215,411 ,21X,4F 10.0) 

1 READ(40, 100)  W1 , W2 , I , J , TYPE , 11 , 12 , 13 , 14 , RT , ZT , BC1 , BC2 , BC3 , BC4 
IF(I.EQ.O)  GO  TO  5 

KK-0 

101  IF(J.LE.JMAX)  GO  TO  2 
ERR-  ERR+1 
WRITE (41, 200) I, J 

2 CT-I4+1 0*13+1 00*12+1 000*11 

IF(I1+I2+I3+I4.NE.0) WRITE (41 , 203)1 ,J, II ,BC1 , I2.BC2 , I4.BC4 

IF(ABS(BC1)+ABS(BC2)+ABS(BC3)+ABS(BC4) .EQ.O. ) GO  TO  3 

UF(IC)  - BC1 

WF(IC)  - BC2 

TANF(IC)-BC4 

CT-CT+IC*10000 

IC  - IC+1 

IF(IC.LE.IOO)  GO  TO  3 
ERR  - ERR+1 
WRITE (41, 201) 

GO  TO  19 

3 CT  - CT+1000000*TYPE 

IF( ABS ( RT )+ABS ( ZT) . NE . 0 . ) CT-CT+1 00000000 
Jl-J-JMIN+1 

IMIN(Jl)  - MIN0(IMIN(J1),I) 

IMAX(Jl)  - MAXO ( IMAX ( J 1 ) , I ) 

IF(W1.NE. POLAR)  GO  TO  4 
RAD  - RT 

ANGLE-ZT/57. 2957795 
RT-  RAD*COS( ANGLE) 

ZT-  RAD*SIN( ANGLE) 

4 WRITE(45)  I,J1 ,RT,ZT,CT 
ICODNT-ICOONT+1 
IF(KK.EQ.l)  GO  TO  404 
IF(W2.NE .LINE)  GO  TO  401 

C****  CALCULATION  OF  POINTS  ON  A STRAIGHT  LINE 
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57300 

57400 

57500 

57600 

57700 

57800 

57900 

58000 

58100 

58200 

58300 

58400 

58500 

58600 

58700 

58800 

58900 

59000 

59100 

59200 

59300 

59400 

59500 

59600 

59700 

59800 

59900 

60000 

60100 

60200 

60300 

60400 

60500 

60600 

60700 

60800 

60900 

61000 

61100 

61200 

61300 

61400 

61500 

61600 

61700 

61800 

61900 

62000 

62100 

62200 

62300 

62400 


READ(40,103)  W1,W3,II,JJ,RT1,ZT1 
IF(W1.NE. POLAR)  GO  TO  402 
RAD  - RT1 

ANGLE  - ZT1/57. 2957795 
RT1»  RAD*COS( ANGLE) 

ZT1-RAD*S IN (ANGLE) 

W1  - 0. 

402  NSTEPS-MAXO(IABS(II-I)  .IABS(JJ-J)) 

ISTEP“( II-I ) /NSTEPS 

J STEP* ( J J - J ) /NSTEPS 
DD-NSTEPS 
DR  - (RT1-RT)/DU 
DZ-  (ZT1-ZT)/DD 
DO  404  N-l, NSTEPS 
IJK  - N 

I- I+ISTEP 
J-J+JSTEP 
JW-JMIN+1 

IMAX( J1 )-MAX0 ( IMAX( J1 ) , I ) 

IMIN( J1 )-MIN0( IMIN( J1 ) , I) 

RT-RT+DR 

ZT-ZT+DZ 

IF(W3 .EQ. SAME)  GO  TO  403 

READ(40,102)  I, J, TYPE, XI , 12, 13, I4.BC1 ,BC2,BC3 ,BC4 
KK-1 

GO  TO  101 

403  CT-CT 

IF( J . EQ. JMAX)  CT-MOD(CT , 1000000 )+l 30000000 
WRITE(45)  I,J1 ,RT,ZT,CT 
ICOUNT  -ICOUNT+1 

IF(I1*I2+I3+I4.NE.0)  WRITE(41 ,203)  I , J , II ,BC1 , I2.BC2, I4.BC4 

404  CONTINUE 
401  GO  TO  1 

C***  STORE  COORDINATES  AND  CODE  IN  I,J  ARRAYS 
5 REWIND  45 

51  READ(45)  I , J1 , RT, ZT,CT 
ICOUNT-ICOUNT-1 

II- I-IMIN(J1)+1 
R(U,J1)-RT 
Z(I1,J1)-ZT 

IF  (I  .EQ.  IMAX(Jl))  CT  - MOD (CT, 1000000)  + 130000000 
C0DE(I1,J1)-CT 
IF ( ICOUNT. GT.O)GO  TO  51 
JRAN  - JMAX-JMIN+1 
C***  DETERMINATION  OF  BAND  WIDTH 
MAXBAN  - 80 
J2"  JRAN-1 
BW-  0 

DO  52  J-1.J2 

NEQ(J)-  2*( IMAXC J)+1-IMIN( J ) ) 

NBAND-  2*(IMAX(j)+3-IMIN(J+l)) 

IF (NBAND.LE. MAXBAN)  GO  TO  52 


f* 


l 
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62500 

WRITE (41 ,207 ) J1 .MAXBAN 

62600 

52 

BW-  MAXO ( BW , NBAND ) 

62700 

NEQ(JRAN)-2*(IMAX(JRAN)+1-IMIN(JRAN)) 

62800 

DO  6 J-l.JRAN 

62900 

IF(IMAX(J)-IMIN(J).LE. MAX-1)  GO  TO  6 

63000 

ERR-  ERR+1 

63100 

J1  - J+JMIN-1 

63200 

WRITE (41, 202)  Jl 

63300 

6 

CONTINUE 

63400 

IF(ERR.NE.O)  GO  TO  19 

63500 

C**** 

CALCULATE  COORDINATES  OF  INTERIOR  GRID  POINTS 

63600 

IF( JRAN.LE.2)  GO  TO  9 

63700 

J2-  JRAN-1 

63800 

DO  8 N-1,500 

63900 

IJK  - N 

64000 

RESID-  0. 

64100 

DO  7 J-2.J2 

64200 

IRAN  - IMAX(J)-IMIN(J) 

64300 

IN  - IMIN( J)-IMIN( J+l )+l 

64400 

IM  - IMIN( J)-IMIN( J-l )+l 

64500 

DO  7 1-2, IRAN 

64600 

IN  - IN+1 

64700 

IM  - IM+1 

64800 

IF( CODE ( I, J).GE. 100000000)  GO  TO  7 

64900 

DR  - (R(IN,J+l)+R(IM,J-l)+R(I+l,J)+R(I-l,J))/4.-R(l,J) 

65000 

DZ-  (Z(IM,J-l)+Z(IN,J+l)+Z(I+l,J)+Z(I-l,J))/4.-Z(I,J) 

65100 

R(I, J)  - R(I,J)+1.8*DR 

65200 

Z(I,J)  - Z(I,J)+1.8*DZ 

65300 

RESID  - RESID  + ABS(DR)+ABS(DZ) 

65400 

7 

CONTINUE 

65500 

IF(N.EQ.l)  RES1-  RESID 

65600 

IF(RESID/RESl.LT.l.E-5)  GO  TO  9 

65700 

8 

CONTINUE 

65800 

c**** 

OUTPUT  OF  NODAL  POINT  COORDINATES 

65900 

WRITE (43, 799) 

66000 

799 

FORMAT(6X,°I° ,3X,°J° ,12X,®MAT  NO.® ,12X,°AREA° ,/) 

66100 

9 WRITE (41, 205 )N 

66200 

DO  16  J-l.JRAN 

66300 

I1»IMIN(J)-1 

66400 

Jl-  JMIN-l+J 

66500 

IRAN  - IMAX( J)-I1 

66600 

DO  16  1-1 , IRAN 

66700 

11-  11+1 

66800 

c*** 

FIX  RADIAL  DISPLACEMENT  IF  R EQUALS  ZERO 

66900 

IF(R(I,J).EQ.0.)CODE(I,J)-(CODE(I,J)/10000)*10000+1000+MOD(CODE( 

67000 

3 

l I, J), 1000) 

67100 

c*** 

SET  TYPE  EQUAL  30  IF  I EQUAL  IRAN  OR  J EQUALS  JMAX 

67200 

IF ( I . EQ . IRAN . OR . J . EQ . JRAN ) CODE ( I , J ) -CODE ( I , J ) / 1 00000000* 1 00000000 

67300 

U ♦ 30000000  +MOD ( CODE ( I, J), 1000000) 

67400 

WRITE (41 ,206)  11 , Jl ,R(I , J) ,Z(I, J) ,CODE(I, J) 

67500 

IF(CODE(I,J).EQ.O)  GO  TO  225 

67600 

TYPE-(CODE(I, J)-100000000  )/1000000 
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67700 

67800 

67900 

68000 

68100 

68200 

68300 

68400 

68S00 

68600 

68700 

68800 

68900 

69000 

69100 

69200 

69300 

69400 

69500 

69600 

69700 

69800 

69900 

70000 

70100 

70200 

70300 

70400 

70500 

70600 

70700 

70800 

70900 

71000 

71100 

71200 

71300 

71400 

71500 

71600 

71700 

71800 

71900 

72000 

72100 

72200 

72300 

72400 

72500 

72600 

72700 

72800 


GO  TO  226 

225  TYPE-0 

226  WRITE (46 ) I, J, TYPE 

16  CODE ( I , J ) - MOD ( CODE ( I, J), 100000000) 

1-0 

J-0 

TYPE-0 

WRITE(46)  I, J, TYPE 

C***  CALCULATION  OF  AREAS  OF  TRIANGLES 
DO  18  J-l.JRAN 

11-  IMIN(J) 

12-  IMAX(J) 

Jl-  JMIN-l+J 
DO  18  I-  11,12 

13-  I-Il+l 

IF( CODE ( 13 , J ) / 1000000 . GT . 25 )GO  TO  18 
IN-  I-IMIN( J+l)+l 
RR(1)-  R(I3,J) 

RR(2>-  R(I3+1,J) 

RR(3)-  R(IN+1 , J+l) 

RR(4)-  R(IN,J+1) 

RR ( 5 ) - RR(  1 ) 

ZZ(1)-  Z( 13, J) 

ZZ(2)-  Z(I3+1 , J) 

ZZ(3)-  ZdN+l.J+l) 

ZZ(4)-  Z( IN , J+l ) 

ZZ(5 )-  ZZ(1) 

RK-  (RR(i)+RR(2)+RR(3)+RR(4))/4. 

ZK-  (ZZ(l)+ZZ(2)+ZZ(3)+ZZ(4))/4. 

AREAT-0 . 

ARM-0 ^ 

DO  17  N-1,4 

AREA  - (RR(N+1)-RR(N))*(ZK-ZZ(N) )-(RK-RR(N) )*( ZZ(N+1)-ZZ(N) ) 
RA«(RR(N)+RR(N+l)+RK)/3. 

ARMaRA*AREA+ARM 
IF(AREA.GT.O. ) GO  TO  17 
ERR-  ERR+l 

WRITE (41 ,220)  N.I.Jl 
17  AREAT-AREAT+AREA 

NTT- (CODE (I+1-IMIN(J) ,J)/1000000) 

IF(NTT.NE.18)  GO  TO  18 
VOL-3 . 1415926536*ARM 
WRITE(43 ,999 ) I, Jl, VOL, NTT 
999  FORMAT( 2X , 215 , 5X , # VOL-* , E10 . 4 , II 5 ) 

18  CONTINUE 

IF(TITLE(13) . NE. WORD 1) WRITE (46)  IMAX.IMIN.JRAN 

19  RETURN 

100  FORMAT( 2A1 , 13 , 215 , 411 , IX , 6F10 . 5 ) 

102  F0RMAT(3I5,4I1 ,21X,4F10.5) 

103  FORMAT(2A1,I3,I5,10X,2F10.5) 

200  FORMAT(26HOJ  EXCEEDS  JMAX  ON  CARD  I-I5,4H  J-I5) 

201  FORMAT(53HOMORE  THAN  99  NODES  HAVE  NON  ZERO  BOUNDARY  CONDITIONS) 


72900 

73000 

73100 

73200 

73300 

73400 

73500 

73600 

73700 

73800 

73900 

74000 

74100 

74200 

74300 

74400 

74500 

74600 

74700 

74800 

74900 

75000 

75100 

75200 

75300 

75400 

75500 

75600 

75700 

75800 

75900 

76000 

76100 

76200 

76300 

76400 

76500 

76600 

76700 

76800 

76900 

77000 

77100 

77200 

77300 

77400 

77500 

77600 

77700 

77800 

77900 

78000 


202  FORMAT ( 36HOMORE  THAN  30  NODAL  POINTS  ON  ROW  J-I5  ) 

203  FORMAT(2I5,3(I8, 1PE17 .7)) 

204  FORMAT( 25HOBOUNDARY  CONDITION  ARRAY/ 1 OHO  NODAL  PT10X6HRADIAL19X5HA 

1XIAL  15X7HSLIDING  / 1H  3X 1 H I4X 1 H J 5X4HC ODE 7 X5H VALUE 9 X 

24HCODE7X5HVALUE9X4HCODE7X5HVALUE) 

205  FORMAT ( 30H 1 COORD INATE S CALCULATED  AFTER  I3.11H  ITERATIONS  / 
14X1HI4X1HJ10X1HR14X1HZ14X4HC0DE  /) 

206  FORMAT(2I5,2F15.4,I15) 

207  FORMAT(21H  BAND  WIDTH  OF  ROW  J-I3,  8HEXCEEDS  15) 

220  FORMATC49H  ZERO  OR  NEGATIVE  AREA  IN  TRIANGULAR  ELEMENT  NO.  II, 

126H  OF  QUADRILATERAL  ELEMENT  12, 2H,  12  ) 

END 

SUBROUTINE  TRIAN(A,B,NEQJ,NEQJ1 ,BW) 

C****  this  SUBROUTINE  TRIANGLIZES  A BLOCK  OF  BANDED  EQUATIONS 
IMPLICIT  REAL*8(A-H,0-Z) 

REAL*8  TITLE , WORD 1 , WORD2 , WORD3 

INTEGER  BW.BWl 

DIMENSION  A(120,80),B(120) 

DOUBLE  PRECISION  RATIO 
N3-  NEQJ+NEQJ1 
DO  1 N-1,NEQJ 
Nl-  N+l 

N2-  MIN0(BW+N-1,N3) 

IF(A(N,1) .EQ.O. ) GO  TO  2 
DO  1 I-N1.N2 

11-  I+l-N 

RATIO-  A(N,I1)/A(N,1) 

IF( RATIO. EQ.O.)GO  TO  1 
B(I)-  B(I)-RATIO*B(N) 

Jl-  BW-I1+1 
DO  3 J-l.Jl 

12-  I-N+J 

3 A(I, J)-  A(I,J)-A(N,I2)*RATIO 

1 CONTINUE 
RETURN 

2 WRITE(41 ,201 ) 

201  FORMAT(49H  ZERO  TERM  ON  MAJOR  DIAGONAL  EXECUTION  TERMINATED) 

STOP 

END 

SUBROUTINE  BACSUB(A,B,NEQ,BW, JRAN) 

IMPLICIT  REAL*8(A-H,0-Z) 

REAL*8  TITLE .WORDl .WORD2 , WORD3 
INTEGER  BW 

DIMENSION  A( 1 20 , 80 ) , B( 1 20 ) , NEQ( 50 ) , DISP ( 1 20 ) 

DOUBLE  PRECISION  DP 
DATA  JR1/121/ 

DO  5 J-l.JRAN 
Jl-  JRAN+l-J 
NEQJ-  NEQ(Jl) 

NEQJi-  NEQ( Jl+1 ) 

IF( Jl .EQ. JRAN)  NEQJI-  0 
Nl-  NEQJ+NEQJl 


I'  H 


f 


f 


78100 

78200 

78300 

78400 

78500 

78600 

78700 

78800 

78900 

79000 

79100 

79200 

79300 

79400 

79500 

79600 

79700 

79800 

79900 

80000 

80100 

80200 

80300 

80400 

80500 

80600 

80700 

80800 

80900 

81000 

81100 

81200 

81300 

81400 

81500 

81600 

81700 

81800 

81900 

82000 

82100 

82200 

82300 

82400 

82500 

82600 

82700 

82800 

82900 

83000 

83100 

83200 


C 


— 


1 

2 


3 


4 

5 


1 


2 


3 


_ X 


BACKSPACE  45 

READ(45)  ((A(I1,I2),I2-1,BW),B(I1),I1-1,NEQJ) 

BACKSPACE  45 

DO  2 N-l.NEQJ 

N2-  NEQJ+l-N 

DP-  B(N2) 

DO  1 K-2.BW 
Kl-  N2-1+K 
IF(Kl.GT.Nl)  GO  TO  2 
DP-  DP-A(N2,K)*DISP(K1 ) 

DISP(N2)  - DP/A(N2, 1) 

JR  - JR1-J 
DO  3 K-l.NEQJ 
A(JR,K)-  DISP(K) 

IF(J.EQ. JRAN)GO  TO  5 
DO  4 K-l , NEQJ 
Kl-  NEQ(J1-1)+K 
DISP(Kl)  - A( JR ,K) 

CONTINUE 

RETURN 

END 

SUBROUTINE  INVRT ( A , ACT , DIM) 

INVERSION  OF  SYMMETRIC  MATRIX 
IMPLICIT  REAL*8 ( A-H , O-Z ) 

REAL*8  TITLE , WORD1 ,WORD2 ,WORD3 

INTEGER  ACT, DIM 

DIMENSION  A(DIM,DIM) ,LOC(61 ) 

DOUBLE  PRECISION  DP 
ABStt)-DABS(X) 

DO  1 N-l.ACT 
LOC(N)-N 
DO  6 Nl-l.ACT 
M-0 

PIVOT-O. 

DO  2 N2-N1 , ACT 
NN-LOC(N2) 

IF  (ABS(A(NN,NN)) .LE.ABS(PIVOT))  GO  TO  2 
M-N2 

PIVOT-A( NN , NN ) 

CONTINUE 

IF  (M.EQ.O)  GO  TO  8 
N-LOC(M) 

LOC(M)-LOC(Nl) 

L0C(N1)-N 
A(N,N)— 1. 

DO  3 J-l.ACT 
A( N , J )-A( N , J ) /PIVOT 
DO  5 Il-l.ACT 
I-LOC(Il) 

DP-A(I.N) 

IF  (N.EQ.I.OR.A(I.N).EQ.O.)  GO  TO  5 
DO  4 Jl-Il , ACT 


.... 


83300 

J-LOC(Jl) 

r 

83400 

IF  (N.EQ.J)  GO  TO  4 

ll 

83500 

A(I,J)-A(I,J>-  A(N, J)*DP 

83600 

A(J,I)-A(I,J) 

r » 

83700 

4 CONTINUE 

83800 

5 CONTINUE 

Li 

83900 

DO  6 I-l.ACT 

84000 

6 A(I,N)«A(N,I) 

|j 

84100 

DO  7 I-l.ACT 

u 

84200 

DO  7 J-l.ACT 

84300 

7 A(I,J)—  A(I,J) 

i : 

84400 

RETURN 

Ll 

84500 

8 WRITE(41 ,9) 

84600 

9 FORMAT  (42H0MATRIX  IS  SINGULAR  - EXECUTION  TERMINATED  ) 

I i 

84700 

CALL  EXIT 

84800 

RETURN 

LJ 

84900 

END 

85000 

SUBROUTINE  CLEAN 

85100 

RETURN 

Ll 

85200 

END 

85300 

SUBROUTINE  TEMPT (TOPT) 

85400 

IMPLICIT  REAL*8(A-H,CKZ) 

85500 

REAL*4  RI , ZI , TIMP , TXME , TXEMP 

85600 

REAL*8  TITLE, WORDl .WORD2 , WORD3  .ISTUFF 

85700 

INTEGER  TOPT, CODE, ERR 

85800 

LOGICAL  L1.L2 

85900 

COMMON  ISTUFF , JMIN , JMAX , ERR , MAX , R , Z , CODE , TITLE , IMAX , IMIN 

86000 

DIMENS ION  RI ( 2000 ) , ZI ( 2000 ) , TIMP( 2000 ) 

86100 

DIMENSION  TEMP( 2000 ) ,RT( 2000) ,ZT(2000) ,D(5) , RR(5) ,ZZ(5),T(5) 

86200 

1 ,ISTUFF(300) ,R(30,50) ,Z(30,50) ,CODE(30,50) ,TITLE(13) 

86300 

2 ,IMAX(50) ,IMIN(50) ,DT2(30) 

86400 

3 ,TTEMP(30,50),IPRINT(30) 

86500 

COMMON/TEM/A( 120,80) 

86600 

EQUIVALENCE  (A(3001 ) ,RT) ,(A(5001 ) ,ZT) , (A(7001 ) .TEMP.TTEMP) 

86700 

DIMENSION  TXEMP(30,50) 

86800 

SIN(X)-DSIN(X) 

86900 

COS(X)-DCOS(X) 

87000 

ABS(X)-DABS(X) 

87100 

ATAN(X)-DATAN(X) 

87200 

READ (40, 100)  TO 

87300 

WRITE (41, 201 )T0 

[j 

87400 

ITOPT-TOPT+1 

Ll 

87500 

GO  T0(20,l ,3,4,41 ) , ITOPT 

87600 

1 READ(40, 100)  T1 

r 

87700 

DT-  T1-T0 

87800 

WRITE (41, 202 )DT 

l- « 

87900 

JRAN1-  JMAX- JMIN 

r 

88000 

DO  2 J-l , JRAN1 

88100 

11-  IMIN(J) 

L. 

88200 

I2-IMAX( J)-l 

88300 

DO  2 1-11,12 

r 

88400 

13-  I+l-Il 

\ . 


88500 

88600 

88700 

88800 

88900 

89000 

89100 

89200 

89300 

89400 

89500 

89600 

89700 

89800 

89900 

90000 

90100 

90200 

90300 

90400 

90500 

90600 

90700 

90800 

90900 

91000 

91100 

91200 

91300 

91400 

91500 

91600 

91700 

91800 

91900 

92000 

92100 

92200 

92300 

92400 

92500 

92600 

92700 

92800 

92900 

93000 

93100 

93200 

93300 

93400 

93500 

93600 


IF( CODE ( 13 , J)/ 1000000. GT. 25)  GO  TO  2 
WRITE (42)  DT 

2 CONTINUE 
RETURN 

3 READ(40, 101 ) NOTEMP, Z0,D4 

READ(42)  (RI(N),ZI(N),TIMP(N), N-l, NOTEMP) 

DO  999  N-l, NOTEMP 
RT(N)-RI(N) 

ZT(N)-ZI(N) 

999  TEMP(N)«TIMP(N) 

GO  TO  5 

4 READC40 , 101 ) NOTEMP, Z0.D4 

READC40, 102)  (RT(N) ,ZT(N) ,TEMP(N) , N-l .NOTEMP) 

GO  TO  5 

41  READ(40, 101 )NDIST 
DO  42  Nl-1 , NDIST 
READ(42)  TXME.TXEMP 
TIME-TXME 
DO  442L3-1 , 30 
DO  442  L4-1.50 

442  TTEMP(L3,L4)-  TXEMP(L3,L4) 

42  CONTINUE 
GO  TO  602 

5 DO  6 N-l, NOTEMP 

6 ZT(N)»  ZT(N)-Z0 
602  JRAN1-  JMAX-JMIN 

D4-D4*D4 

IF ( TOPT . NE . 4 ) WRITE ( 4 1 , 203 ) 

IF( TOPT . EQ . 4 ) WRITE (4 1 , 205 )TIME 
DO  16  J-1.JRAN1 

11-  IMIN(J) 

12-  IMAX(J)-1 
Jl-JMIN-l+J 

DO  1601  1-11,12 

13-  I+l-Il 

IF ( CODE ( 13 , J ) / 1000000 . GT . 25 ) GO  TO  1601 
IN-  I-IMIN(J+1)+1 

RK-  (R(I3,J)+R(I3+l,J)+R(IN,J+l)+R(IN+l,J+l))/4. 
ZK-  (Z(I3, J)+Z(I3+1 , J)+Z(IN, J+l )+Z( IN+1 , J+l) )/4. 
IF(ITOPT.NE.5)GO  TO  1602 
RR(1)-R(I3,J) 

RR(2)-R(I3+1,J) 

RR(3)-R(IN+1 , J+l) 

RR(4)-R(IN, J+l) 

ZZ(1)-Z(I3, J) 

ZZ(2)-Z( 13+1 , J) 

ZZ(3)-Z( IN+1, J+l) 

ZZ(4)-Z(IN, J+l) 

T(1)-TTEMP(I3,J) 

T(2)-TTEMP(I3+1 , J) 

T(3)-TTEMP(IN+1,J+1) 

T(4)-TTEMP(IN, J+l ) 


•A 
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E 


5 


93700  GO  TO  1201 

93800  1602  14-13 

93900  15  - IN 

94000  IF  (I  .LT.  (II  ♦ I 2 ) / 2 ) GO  TO  601 

94100  14  - 14  ♦ 1 

94200  15  - 15  + 1 

94300  601  THETA-1.570795 

94400  IF(ABS(R(I5,J+l)-R(I4,J)).GT..001)THETA-ATAN((Z(I5,J+l)-Z(I4,J))/( 

94500  1 R(I5,J+1)-R(I4,J))) 

94600  D(l)-D4 

94700  D(2)-D4 

94800  D(3)-D4 

94900  D(4)-D4 

95000  DO  12  N-l, NOTEMP 

95100  AA-RT(N )-RK 

95200  BB-ZT(N)-ZK 

95300  DD-AA*AA+BB*BB 

95400  IF(DD.GE.D4)G0  TO  12 

95500  LI -.TRUE. 

95600  L2- . TRUE . 

95700  IF( AA*SIN ( THETA ) -BB*COS ( THETA) . LT . 0 . ) L 1 - . FALSE . 

95800  IF(AA*COS(THETA)*BB*SIN(THETA) .LT.O. )L2-. FALSE. 

95900  DO  11  L-1,4 

96000  GO  TO(7,8,9,10),L 

96100  7 IF( .NOT. LI .OR. .NOT.L2)  GO  TO  11 

96200  GO  TO  10 

96300  8 IF(L1 .OR. .NOT.L2)  GO  TO  11 

96400  GO  TO  10 

96500  9 IF(L1.0R.L2)  GO  TO  11 

96600  GO  TO  10 

96700  10  IF(DD.GE.D(L))GO  TO  12 

96800  D(L)-DD 

96900  T(L)-TEMP(N) 

97000  RR(L)-  RT(N) 

97100  ZZ(L)-  ZT(N) 

97200  GO  TO  12 

97300  11  CONTINUE 

97400  12  CONTINUE 

97500  1201  T( 5 )— T( 1 ) 

97600  CC-  0. 

97700  DTI-  0. 

97800  RR(5 )-  RR(1) 

97900  ZZ(5)-  ZZ(1) 

98000  D(5)-D(l) 

98100  DO  17  N-l, 4 

98200  IF(ITOPT.EQ.5)GO  TO  1701 

98300  IF(D(N) .GE.D4.0R.D(N+1) .CE.D4)  GO  TO  17 

98400  1701  AJ-  RR(N+1 )-RR(N) 

98500  BJ-  ZZ(N*1)-ZZ(N) 

98600  AK-  RK-RR(N) 

98700  BE-  ZK-ZZ(N) 

98800  A1EA-AJ*BK-AK*BJ 
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98900 

99000 

99100 

99200 

99300 

00100 

00200 


00100 

00200 

00100 

00200 

00300 

00400 


00100 


00100 

00200 


00100 

00200 

00300 

00400 

00500 

00600 

00700 

00800 

00900 

01000 

00100 

00200 

00300 

00400 

00300 

00600 

00700 

00800 

00900 


WRITE(41 ,599)  AJ,BK, AK,BJ, AREA 
599  FORMAT(2X,°AJ,BK,AK,BJ,AREA-#,5E13.6) 

IF(AREA.EQ.O. ) GO  TO  17 
C-  ZZ(N+1 )-ZK 
DX-  RK-RR(N+1 ) 

COMM-  (RR(N)+RR(N+l)+RK)/6./AREA 

DTI-  C0MM*(BK*BJ+AJ*AK)*T(N+1)+C0MM*(C*BJ-DX*AJ)*T(N)+DT1 


CC-CC+CO*M*( BJ*BJ+AJ*AJ ) 
17  CONTINUE 

DTI-  DT1/CC-T0 
WRITE (42)  DTI 
IPRINT(I3)-I 
1601  DT2(I3)-DT1 


IRAN-12+1-11 


16  WRITE (41 ,204)  J1 , ( IPRINT(I) ,DT2(I) , 1-1 , IRAN) 
20  RETURN 


100  FORMAT(F10.5) 

101  FORMATdlO , 2F10.5) 

102  FORMAT(2F10.5,E15.6) 

201  TORMAT( 1H1 , 20X.48HT  EMPERATURE  DISTRIBUTION 

1 / / 10X23HREFRENCE  TEMPERATURE  IS  F10.3) 

202  FORMAT( 6SH0THE  TEMPERATURE  DIFFERENCE  IS  UNIFORM  THROUGHOUT  THE  BO 
IDT  AMD  IS  F10.3  ) 

203  FORMAT ( 64H0THE  TEMPERATURE  DIFFERENCES  WERE  DETERMINED  FROM  AN  INP 
1UT  TABU) 

204  FORMAT( 7HOROW  J-I3/(3H  I-I3.2X3HDT-F8.2, 

1 4X2HI-I3 , 2X3HDT-F8 . 2 , 

1 4X2HZ-I3 , 2X3HDT-F8 . 2 , 

1 4X2HI-I3 , 2X3HDT-F8 . 2 , 

1 4X2HI-I3 , 2X3HDT-F8 . 2 ) ) 

205  FORMAT ( 7 2H0THE  TEMPERATURE  DIFFERENCES  WERE  DETERMINED  FROM  THE  DI 
1 STRIBUTION  AT  T-,1PE12.5,19H  AS  FOUND  BY  AMG065) 

END 

SUBROUTINE  STRESS(U,W,R,Z,I,J,I3,J3,CODE,IC,IN,SSMAX,SSMIN,IJSS) 
C***«*  CALCULATION  OF  STRESSES  AND  STRAINS  IN  A QUADRILATERAL 
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01000 

01100 

01200 

01300 

01400 

01500 

01600 

01700 

01800 

01900 

02000 

02100 

02200 

02300 

02400 

02500 

02600 

02700 

02800 

02900 

03000 

03100 

03200 

03300 

03400 

03500 

03600 

03700 

03800 

03900 

04000 

04100 

04200 

04300 

04400 

04500 


04600 

04700 

04800 

04900 

05000 

05100 

05200 

05300 

05400 

05500 

05600 

05700 

05800 

05900 

06000 

06100 


C DISPLACEMENT  FIELD  IS  FITTED  BY  A FOUR  TERM  PARABALOID 

C TO  THE  FOUR  CORNER  AND  CENTER  POINT  DISPLACEMENTS  BY 

C LEAST  SQUARES 

IMPLICIT  REAL*8(A-H,0-Z) 

REAL*8  TITLE , WORD1 , WORD2 , WORD3 
DOUBLE  PRECISION  F2,S28 
INTEGER  CODE 
LOGICAL  EL1.EL2 

EQUIVALENCE  (EP.EPR) ,(EP(2) ,EPT) , (EP(3) ,EPZ) , (EP(4) .GAMRZ) ,(SIG, 
1SIGR) , ( SIG(2) , SIGT) , (SIG(3) , SIGZ) , (SIG(4) ,TAURZ) 

EQUIVALENCE  (SS ,SIG) , (SS(5 ) , SIGMAX) , (SS(6) , SIGMIN) , (SS(7) .TAUMAX) , 
1 ( SS(8 ) ,EP) , ( SS( 1 2 ) , EPMAX) , ( SS( 13) ,EPMIN) , (SS( 14) .GAMMAX) 

DIMENSION  PHI(4,5),PTP(4,4),PU(4),PW(4),F2(2),S28(2,8) 

1 ,U1(5 ) ,W1(5) , S4(4,4) ,ALF(4) ,BET(4) ,U(30, 50) ,W(30,50) 

2 ,C(4 ,4) ,EP(4) ,ET(4) , SIG(4) .CODE (30, 50) , R(30,50) ,Z(30, 50) 

3 , SSMAX( 14 ) , SSMIN ( 14 ) , SS ( 14 ) , I JSS( 14 , 4 ) 

DATA'  PHI( 1,1) ,PHI(1 ,2) , PHI( 1 ,3) , PHI( 1 ,4) ,PHI( 1,5)  /5*1./, 

1PHI(2 , 1 ) ,PHI(3 , 1)  , PHI (4, 0/3*0. 0 / 

COS(X)-DCOS(X) 

SIN(X)-DSIN(X) 

ABS(X)-DABS(X) 

SQRT(X)-DSQRT(X) 

ATAN(X)-DATAN(X) 

SIGN(X,Y)-DSIGN(X,Y) 

RE AD (44)  C , ET , RR , ZZ , F2 , S28 
U1(2)-U(I,J) 

U1(3)-U(I+1,J) 

U1(4)-U(IN,J+1) 

U1(5)-U(IN+1,J+1) 

W1(2)-W(I,J) 

W1(3)-W(I+1,J) 

W1(4)-W(IN,J+1) 

W1(5)*W(IN+1 ,J+1) 

PHI(2,2)-R(I,J)-RR 
PHI( 2 , 3 ) *R( I+l , J )-RR 
PHI(2,4)-R(IN,J+1)-RR 
PHI(2,5)-R(IN+1,J+1)-RR 
PHI(3,2)-Z(I,J)-ZZ 
PHI(3 , 3)*Z( I+l , J)-ZZ 
PHI(3,4)*Z(IN, J+1)-ZZ 
PHI(3,5)-Z(IN+1,J+1)-ZZ 
DO  1 M-2,5 

1 PHI(4,M)-PHI(2,M)**2  ♦ PHI(3,M)**2 
DO  2 M-1,2 

DO  2 MM-2,5 

2 F2(M)-F2(M)-S28(M, 2*MM-3)*U1 (ltO-S28(M, 2*MM-2)*W1 (MM) 

Ul(l)-F2(l) 

Wl(l)-F2(2) 

DO  4 N-1,4 
DO  3 M-1,4 
S4(M,M)-0. 

DO  3 MM-1,5 
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1. 


0 

0 

L 
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1 

I 

1 
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I 

[ 
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06200 

06300 

06400 

06500 

06600 

06700 

06800 

06900 

07000 

07100 

07200 

07300 

07400 

07500 

07600 

07700 

07800 

07900 

08000 

08100 

08200 

08300 

08400 

08500 

08600 

08700 

08800 

08900 

09000 

09100 

09200 

09300 

09400 

09500 

09600 

09700 

09800 

09900 

10000 

10100 

10200 

10300 

10400 

10500 

10600 

10700 

10800 

10900 

11000 

11100 

11200 

11300 


3 S4(N,M)-S4(N,M)+PHI(N,MM)*PHI(M,MM) 

PU(N)-0. 

PW(N)-0. 

DO  4 M-1,5 

PU(N)-PU(N)+PHI(N,M)*U1(M) 

4 PW(N)-PW(N)+PHI(N,M)*W1(M) 

CALL  INVRT(S4,4,4) 

DO  5 N-1,4 
ALF(N)-0. 

BET(N)-0. 

DO  5 M-1,4 

ALF(N)-ALF(N)+S4(N,M)*PU(M) 

5 BET(N)-BET(N)+S4(N,M)*PW(M) 

EPR-ALF(2) 

EPT  • U1(1)/RR 
EPZ-BET(3) 

GAMRZ-ALF(3)+BET(2) 

ANGLE" . 5*ATAN ( GAMR2/ ( EPR-  EPZ))*57.2958 

IF  (EPR  .LT.  EPZ)  ANGLE  - ANGLE  ♦ SIGN(90 .ODO.GAMRZ) 
TEM-(EPR  +EPZ  )/2. 

TEM1  -SQRT(((  EPR-  EPZ)/2 . )**2+GAMRZ**2/4. ) 

EPMAX-TEM+TEMl 
EPMIN-TEM-TEM1 
GAMMAX-2 . *TEM1 
DO  8 N-1,4 
SIG(N)— ET(N) 

DO  8 M-1,4 

8 SIG(N)-  SIG(N)+C(N,M)*EP(M) 

TEM  - (SIGR+SIGZ)/2 . 

TEM1  - SQRT( ( ( SIGR-SIGZ) /2 . )**2+TAURZ**2 ) 

SIGMAX-  TEM+TEM1 
SIGMIN-  TEM-TEM1 
TAUMAX  - TEM1 
IF  (MOD(IC, 19))  7,6,7 

6 WRITE(41 , 100) 

7 WRITE(41, 101)13, J3,RR,ZZ,SIGR,SIGT,SIGZ,TAURZ, SIGMAX, SIGMIN, 
1TAUMAX , ANGLE , EPR , EPT , EPZ , GAMRZ , EPMAX , EPMIN , GAMMAX , CODE ( I , J ) 

IF(CODE(I,J).NE. 18000000)  GO  TO  661 

WRITE(47 ,499)  13, J3,SIGR,SIGT, SIGMAX, SIGMIN, CODE ( I, J) 

499  FORMAT(1HOI3,I4,1P4E13.4,I15) 

661  IC-IC+1 

DO  11  K-1,14 

IF(SS(K) . GT. SSMIN(K) )GO  TO  10 
SSMIN(K)-SS(K) 

IJSS(K,1)-I3 

IJSS(K,2)-J3 

10  IF(SS(K).LT.SSMAX(K))GO  TO  11 
SSMAX(K)-SS(K) 

IJSS(K,3)-I3 

IJ88(K,4)-J3 

11  CONTINUE 
ELI*. FALSE. 


11 

li 


1. 


11400 

EL 2-. FALSE. 

11500 

IF  (CODE(IN, J+l )/1000000  .LT.  26)  ELI  - .TRUE. 

11600 

IF  (C0DE(I+1,J  )/1000000  .LT.  26)  EL2  - .TRUE. 

11700 

IF  ( .NOT. ELI .AND.. NOT. EL2)  N-l 

11800 

IF  (.NOT. ELI. AND. EL2)  N - 2 

11900 

IF  (ELI .AND.. NOT. EL2)  N-3 

12000 

IF  (ELI .AND.EL2)  N - 4 

12100 

WRITE (46 ) I , J , RR , ZZ , SIGR, SIGZ , SIGT , TAURZ , SIGMAX , SIGMIN, TAUMAX 

12200 

1 , EPR , EPZ , EPT , GAMRZ , EPMAX , EPMIN , GAMMAX , N 

12300 

100  FORMAT ( 9H1  . J/ 

12400 

1 5X  1 1 HCOORDINATE S37X33HS  TRESSES 

/ S T R 

12500 

1A  I N S /8H  ANGLE7X1HR8X1HZ4X9HRADIAL  R2X10HHOOP  THETA5X8HAXIA 

12600 

2L  Z3X10HSHEAR  R-Z6X7HMAXIMUM6X7HMINIMUM4X9HMAX  SHEAR  ) 

12700 

101  FORMAT(1HOI3,I4,OPF8.3,F9.3,1P7E13.4/1H  0PF7.2.17X1P7E13 

.4,115) 

12800 

102  FORMAT  (1H  1PUE11.3) 

12900 

RETURN 

13000 

END 

13100 

SUBROUTINE  PRESBC(ICF,JCF,H,TE, SIDE, LC, ERR) 

13200 

C***  THIS  SUBROUTINE  READS  THE  PRESSURE  BOUNDARY  CONDITION 

DATA 

13300 

C DATA  MAY  BE  READ  IN  FOR  AN  ELEMENT  SIDE  AT  A TIME  OR  FOR  A LINE 

13400 

C WHICH  HAS  I OR  J AS  A CONSTANT.  THE  SUBROUTINE  MAY  BE 

CHANGED 

13500 

C TO  FIT  PARTICULAR  NEEDS. 

13600 

IMPLICIT  REAL*8 ( A-H , O-Z ) 

13700 

REAL*8  TITLE .WORD1 , WORD2 , WORD3 

13800 

INTEGER  SIDE,  SIDET,  ERR 

13900 

C 

14000 

DIMENSION  ICF(200),  JCF(200),  H(200), 

14100 

1 TE ( 200 ) , SIDE (200) 

14200 

C 

14300 

LC-0 

14400 

1 READ(40, 100)  11 , J1 , 12, J2,HT,TET, SIDET 

14500 

IF( I 1.EQ.0) RETURN 

14600 

LC-LC+1 

14700 

IF(LC.GT. 200)  GO  TO  5 

14800 

ICF(LC)-I1 

14900 

JCF(LC)-J1 

15000 

H(LC)-HT 

15100 

TE(LC)-TET 

15200 

SIDE (LC) -SIDET 

15300 

IF(I2.EQ.O)  GO  TO  1 

15400 

NSTEPS-MAXO( IABS( 12-11 ) , IABS( J2-J1 ) ) 

15500 

ISTEP-( I2-I 1 ) /NSTEPS 

15600 

JSTEP-( J2-J 1 ) /NSTEPS 

15700 

DO  3 N-l, NSTEPS 

15800 

LC-LC+1 

15900 

IF(LC.GT. 200)  GO  TO  5 

16000 

ICF(LC)-ICF(LC-1 )+ISTEP 

16100 

JCF(LC)-JCF(LC-1 )+JSTEP 

16200 

H(LC)-HT 

16300 

TE(LC)-TET 

16400 

3 SIDE (LC) -SIDET 

16500 

GO  TO  1 

I 

I 
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Hi 


16600 

5 

ERR-ERR+1 

16700 

WRITE (41, 209) 

16800 

RETURN 

16900 

100 

FORMAT(4I5 , 2F10 . 5 , 15 ) 

*7000 

209 

FORMAT ( 7 OHONUMBER  OF  ELEMENT  SIDES  WITH  PRESSURE  BOUNDARY  CONDITIO 

17100 

INS  EXCEEDS  200) 

17200 

END 

17300 

SUBROUTINE  SETUP(A.B) 

17400 

c*** 

ASSEMBLE  STIFFNESS  MATRIX  OF  STRUCTURE  IN  THE  FORM  OF  A BAND 

17500 

c 

EQUATIONS  ARE  MODIFIED  FOR  BOUNDARY  CONDITIONS  AND 

17600 

c 

TRIANGLIZED  BEFORE  BEING  WRITTEN  ON  TAPE 

17700 

IMPLICIT  REAL*8(A-H,0-Z) 

17800 

REAL*8  TITLE , WORD1 , WORD2 , WORD3 

17900 

INTEGER  BW, BCCODE , SLCODE , UCODE , WCODE , CODE , BW1 , ERR, PCODE 

18000 

COMMON  BC,  JMIN.JMAX, ERR, MAX, R,Z, CODE, TITLE, IMAX.IMIN, 

18100 

1 CONPR, IP, JP.P.TAU, PCODE, BW.NEQ.JRAN,  S,F,NPCARD 

18200 

2 , ITOPT , NP , MN 

18300 

EQUIVALENCE(IFIX(1 ) , UCODE) ,( IFIX(2) , WCODE) ,(BC( 1 ) ,UF) , (BC( 101 ) ,WF) 

18400 

1 , (BC(201 ) ,TANF) 

18500 

DIMENSION  R(30,50) ,Z(30,50) , CODE (30,50 ) ,TITLE(1 3) ,UF(100),WF(100) 

18600 

1 , TANF( 100) ,IMAX(50) ,IMIN(50) ,BC( 100,3) ,CONPR(16, 15) ,IP(200) 

18700 

2 , JP(200) ,P(200) ,TAU(200) ,PCODE(200) ,NEQ(50) , S(8,8) 

18800 

3 ,A( 120,80) ,B( 120) , IFIX(2) ,F(8) 

18900 

NP-1 

19000 

REWIND  45 

19100 

IPT-0 

19200 

DO  1 1-1,120 

19300 

B(I)-  0. 

19400 

DO  1 J-1,80 

19500 

1 

A(I, J)-  0. 

19600 

DO  20  J-l.JRAN 

19700 

IRAN-  NEQ(J)/2 

19800 

IF(J.EQ.JRAN)  GO  TO  6 

19900 

DO  501  I-l.IRAN 

20000 

-.1-  CODE  ( I , J ) / 1000000 

20100 

IF( MN1 . GT . 0 . AND . MN1 . LT . 26 )MN-MN1 

20200 

IF(MN1 .GT.25 ) GO  TO  501 

20300 

c**+**  CHECK  FOR  PRESSURE  ON  THE  ELEMENT 

20400 

11-  IMI.  (J)t-I-l 

20500 

Jl-  JMIN+J-1 

20600 

IF(NPCARD.LT.NP)GO  TO  1002 

20700 

IS-IP(NP) 

20800 

JS-JP(NP) 

20900 

IFdS.NE.il  .OR.  JS.N_.J1)G0  TO  1002 

21000 

PI-P(NP) 

21100 

PT-TAU(NP) 

21200 

IPT-PCODE (NP ) 

21300 

NP-NP+1 

21400 

1002  CALL  STIFFQ( I , J , PI , PT , IPT ) 

21500 

c*** 

ADD  ELEMENT  STIFFNESS  TO  TOTAL  STIFFNESS 

21600 

Nl-  2MI-1) 

21700 

N2-  2*( I+IMIN( J)-1-IMIN( J+l ) )+NEQ(J) 

i 

H 
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21800 

21900 

22000 

22100 

22200 

22300 

22400 

22500 

22600 

22700 

22800 

22900 

23000 

23100 

23200 

23300 

23400 

23500 

23600 

23700 

23800 

23900 

24000 

24100 

24200 

24300 

24400 

24500 

24600 

24700 

24800 

24900 

25000 

25100 

25200 

25300 

25400 

25500 

25600 

25700 

25800 

25900 

26000 

26100 

26200 

26300 

26400 

26500 

26600 

26700 

26800 

26900 


DO  5 N-1,4 

N3  - N2-N1+1-N 
NN1-  Nl+N 
NN2-  N2+N 

B(NN1)-  B(NN1)+F(N) 

B(NN2)-  B(NN2)+F(N+4) 

LX-  0 

DO  4 L-  H,4 
LX-  LL+1 

A(NN1,LL)-  A(NN1,LL)+S(N,L) 

4 A(NN2,LL)-  A(NN2,LL)+S(N+4,L+4) 

DO  5 L-1,4 

LI-  N3+L 

5 A(NN1,L1)-  A(NNl,Ll)+S(N,L+4) 

501  CONTINUE 

C***  MODIFY  THE  BLOCK  OF  EQUAT  ..S  FOR  FIXITY  AND  APPLIED  LOADS 

6 DO  11  1-1 , IRAN 

BCCODE-  MOD(CODE( I, J), 1000000) 

IF( BCCODE . EQ . 0 ) GO  TO  11 
SLCODE-  MOD(BCC  DE.10) 

WCODE-  MOD( BCCODE, 1000)/ 100 
UCODE-  MOD ( BCCODE, 1 0000 )/1000 
NCODE-  MOD ( BCCODE, 1 000000) /10000 
IF(NCODE.EQ.O)  NCODE-lOO 
IU-  2*(I-l)-*-i 
IW-  IU+1 

IF(UCODE .EQ.2)  B(IU)-  B(IU)+UF(NC0DE)/6. 2831853 
IF(WCODE.EQ.2)  B(IW)-  B( IW) +WF( NCODE )/6. 283 1853 
IF( SLCODE. EQ.O)  GO  TO  8 
ALF-  TANF( NCODE) 

B(IU)«  B(IU)+ALF*B(IW) 

B( IW)-  0. 

A(IU,1)-  «( IU, 1)+ALF*( ALF*(A( IW, 1)+1 . )+2.*A(IU,2)) 

A(IU,2)-  -ALF 
A(IW,1)-  1. 

BW1-  BW-1 . 

DO  7 N-2.BW1 

A( IU,N+1)-  A(IU,N+1)+ALF*A(IW,N) 

A(IW,N)-  0. 

II-  IU+l-N 

IF(II.LT.l)  GO  TO  7 

A(II,N)-  A( II ,N)+ALF*A( II , N+l ) 

A(II,N+1)-  . 

7 CONTINUE 
A(IW,BW)-  0. 

8 NMAX-  NEQ( J )+NEQ( J+l ) 

IF  (J  .EQ.  JRAN)  NMAX  - NEQ(J) 

DO  10  N-1,2 
IR-IU+N-1 

IF(IFIX(N).NE.l)  GO  TO  10 
DO  9 N1-2.BW 
II-  IR+1-N1 
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27000 

IJ-  IR+N1-1 

27100 

IFQI.GT.O)B(II)-  B(II)-A(II,N1)*BC(NC0DE,N) 

27200 

IF( IJ . LE . NMAX)B( IJ)m  B(IJ)-A( IR,N1 )*BC(NCODE,N) 

27300 

IF(II.GT.O)  A(II,Nl)-  0. 

27400 

9 

A(IR,N1)-  0. 

27500 

A(IR,1)-  1. 

27600 

B(IR)«  BC( NCODE, N) 

27700 

10 

CONTINUE 

27800 

11 

CONTINUE 

27900 

IF(J.EQ.JRAN)  GO  TO  17 

28000 

IRANI-  IMAX( J+l )+l-IMIN( J+l ) 

28100 

DO  16  1-1, IRANI 

28200 

BCCODE-  MOD(CODE(I, J+l), 1000000) 

28300 

IF(BCCODE.EQ.O)  GO  TO  16 

28400 

SLCODE-  MOD( BCCODE, 10) 

28500 

WCODE-  MOD(BCCODE,1000)/100 

28600 

UCODE-  MOD(BCCODE, 10000)/1000 

" 700 

NCODE-  MOD (BCCODE, 1000000 )/l0000 

28800 

IF(NCODE.EQ.O)  NCODE-  100 

28900 

IU-  NEQ(J)+2*(I-1)+1 

29000 

IW-  IU+1 

29100 

IF( SLCODE .EQ.O)  GO  TO  13 

29200 

ALF-  TANF( NCODE) 

29300 

BW1-BW-1 

29400 

DO  12  N-2.BW1 

29500 

II-  IU+l-N 

29600 

IF(II .GT.NEQ(J) .OR. II .LT. 1)  GO  TO  12 

29700 

A( II ,N)-  A( II ,N)+ALF*A( II ,N+1 ) 

29800 

A(II,N+1)-  0. 

29900 

12 

CONTINUE 

3C000 

13 

DO  15  N-1,2 

30100 

IF( IFIX(N) . NE . 1 ) GO  TO  15 

200 

IR-IU+N-1 

30300 

DO  14  Nl-2 ,BW 

30400 

II-IR+1-N1 

30500- 

IF(II .LT. 1 .OR. II .GT.NEQ( J ) ) GO  TO  14 

30600 

B(II)-  B(II)-A(II,N1)*BC(NC0DE,N) 

30700 

A( II, Nl)-  0. 

30800 

14 

CONTINUE 

30900 

15 

CONTINUE 

31000 

16 

CONTINUE 

31100 

c*** 

TRIANGLIZE  THE  J BLOCK  OF  COEFFICIENTS 

31200 

17 

NEQJ-  NEQ(J) 

31300 

..KQJ1-  NEQ(J+1) 

31400 

IF( J.EQ. JRAN)NEQJl-0 

31500 

N4  - NEQJ+NEQJ1 

31600 

CALL  TRIAN(A,B,NEQJ,NEQJ1 ,BW) 

31700 

WRITE(45)( (A(N1 ,N2) ,N2-1 ,BW) ,B(N1 ) ,N1-1 ,NEQJ) 

31800 

IF ( J . EQ . JRAN )GO  TO  20 

31900 

DO  18  Nl-l.NEQJl 

32000 

N3-  NEQJ+N1 

32100 

B(N1)-  B(N3) 

i i 


» 


\ 

0 


37400 

37500 

37600 

37700 

37800 

37900 

38000 

38100 

38200 

38300 

38400 

38500 

38600 

38700 

38800 

38900 

39000 

39100 

39200 

39300 

39400 

39500 

39600 

39700 

39800 

39900 

40000 

40100 

40200 

40300 

40400 

40500 


3 S(I,J)-S(I,J)+A(K,I)*CA(K,J) 

DO  4 J-1,4 

4 F(I)«F(I)+A(J,I)*ET(J) 

VOL  - DEL*( RI+RJ+RK) / 6 . 

DO  42  1 - 1,6 

F(I)  ■ F(I)*VOL 
DO  42  J - 1,6 
S(I,J)  - S(l,J)*VOL 
42  S(J,I)  - S(I,J) 

IF  (ABS(P)+ABS(TAU)+ABS(BFR)+ABS(BFZ).EQ.O.)  GO  TO  5 

AJ»  RJ-RI 

AX-  RK-RI 

BJ-  ZJ-ZI 

BK-  ZK-ZI 

IF(ABS(BFR)+ABS(BFZ) .EQ.O. )GO  TO  41 
Xl-(RI+RJ+RK)/6 . 

X2-(RI*( RI+RJ+RK ) +RJ*( RJ+RK )+RK*RK ) / 1 2 . 
X3-((RI*RI+RJ*RJ+RK*RK)*(RI+RJ+RK)+RI*RJ*RK)/20. 

X4-(BK*( RI+2 . *RK+RJ)+BJ*( RI+2 .*RJ+RK) ) /24 . 

X5« ( BK*<  RK*( 2 . * ( RI+R J ) +3 . *RK) +RI*( RI+RJ ) +RJ*RJ ) + 

1 BJ*( RJ*( 2 . *( RI+RK)+3 . *RJ)+RI*( RI+RK)+RK*RK) ) /60 . 

41  F(l)-  F( 1 )+(TAU*AJ-P*BJ)*(RI/2 .+AJ/6 . )+BFR*( (RJ*BK-RK*BJ)*X2+ 

1 (BJ-BK)*X3+(AK-AJ)*X5) 

F ( 2 ) - F ( 2 ) + ( TAU*BJ+P*AJ )*(RI/2.+AJ/6.)+BFZ*(( RJ*BK-RK*BJ ) *X 1 + 

1 ( BJ-BK)*X2+( AK-AJ )*X4 ) 

F ( 3 ) - F( 3 )+( TAU*AJ-P*BJ ) * ( RI / 2 . +AJ / 3 . ) +BFR* ( -RI*BK*X2+BK*X3-AK*X5 ) 
F ( 4 ) - F ( 4 ) + ( TAU*BJ+P*AJ )*( RI/ 2 . +AJ/3 . )+BFZ*( -RI*BK*X 1 +BK*X2-AK*X4 ) 
F( 5 ) - F ( 5 ) +BFR*( RI*B J*X2 -B J*X3+AJ*X5 ) 

F( 6 )»  F(6)+BFZ*( RI*BJ*X1 -3J*X2+AJ*X4 ) 

5 CONTINUE 
RETURN 
END 


